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fH , We review recently developed new powerful techniques to compute a class of Feynman 



diagrams at any loop order, known as sunrise-type diagrams. These sunrise-type topolo- 
gies have many important applications in many different fields of physics and we believe it 
to be timely to discuss their evaluation from a unified point of view. The method is based 
on the analysis of the diagrams directly in configuration space which, in the case of the 
sunrise-type diagrams and diagrams related to them, leads to enormous simplifications as 
compared to the traditional evaluation of loops in momentum space. We present explicit 
formulae for their analytical evaluation for arbitrary mass configurations and arbitrary di- 
mensions at any loop order. We discuss several limiting cases of their kinematical regimes 
which are e.g. relevant for applications in HQET and NRQCD. We completely solve the 
problem of renormalization using simple formulae for the counterterms within dimensional 
regularization. An important application is the computation of the multi-particle phase 
space in D-dimensional space-time which we discuss. We present some examples of their 
numerical evaluation in the general case of L>-dimensional space-time as well as in inte- 
ger dimensions D = Dq for different values of dimensions including the most important 
practical cases Dq = 2, 3, 4. Substantial simplifications occur for odd integer space-time 
dimensions where the final results can be expressed in closed form through elementary 
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functions. We discuss the use of recurrence relations naturally emerging in configuration 
space for the calculation of special series of integrals of the sunrise topology. We finally 
report on results for the computation of an extension of the basic sunrise topology, namely 
the spectacle topology and the topology with an irreducible loop addition. 
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1 Introduction and motivation 



In the present-day analysis of various physical processes the complexity of the problem often 
makes it impossible to find a complete mathematical description of the experimentally measured 
or observed quantities. In order to proceed one usually simplifies the mathematical model which 
describes the bulk of the phenomena. The more delicate effects are accounted for through an 
expansion in terms of a small parameter which has to be defined such that one retains full 
quantitative control in the analysis. The art of a physicist consists in splitting the whole 
analysis into a main part and into small corrections to the main part for which one develops 
a perturbative expansion. One may quote one of the greatest physicists of the last century, 
L.D. Landau, who used to say that physics starts when a small parameter of the problem has 
been found. The usefulness and power of such an approach has been proven over and over for 
many years. An example that comes into mind are the first accurate calculations analyzing 
high precision data from celestial mechanics in terms of Newton's law of gravitation. When 
observational data on the motion of the planets became more accurate it was no longer sufficient 
to take into account only their gravitational interaction with the sun but it was necessary also 
to include the small corrections arising from the mutual gravitational interactions between the 
planets themselves. This was done through perturbation theory since the full equations for 
the planetary movements with mutual interaction fully taken into account are too difficult to 
treat analytically even within standard classical mechanics. A remarkable fact is that these 
corrections are very small but still have to be taken into account since the data are extremely 
precise. This is well illustrated by a historical example taken from astronomy which was the 
most precise scientific disciphne in the past: the planets of the solar system close to Earth were 
directly observed while an eighth planet - Neptune - was first introduced in order to explain 
shght perturbations in the orbit of Uranus. The orbit of Neptune was theoretically calculated 
with high precision and the planet was directly observed later at the predicted place. 

Let us now turn to high energy physics. The main interest of fundamental research in 
high energy physics is to find new physical phenomena and to place them into a consistent 
logical picture of the natural world which mankind inhabits. It is generally believed that the 
structure of fundamental interactions (except gravity) is qualitatively understood at the energies 
available at present accelerators. Basically there are two trends for the search of new physics. 
One is aiming at a direct discovery of new physical phenomena and fundamental constituents of 
matter (elementary particles) by moving to higher and higher energies. This requires ever more 
powerful accelerating facilities while the precision of the experimental detection and theoretical 
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calculations can remain at a moderate level. The other is an indirect search for new physics at 
low energies (mainly through radiative corrections) which is based on very precise experimental 
data at moderate energies and a high accuracy of the theoretical calculations 12 El IH El • 
A very recent example of the second category is the ongoing discussion about the possible 
appearance of deviations from the Standard Model (SM) of strong and electroweak interaction 
in new data on the muon anomalous magnetic moment [HI El El Ej- This two-fold road to 
scientific discovery is well known and has been pursued for ages as illustrated by the above 
example from celestial mechanics. 

The muon anomalous magnetic moment has been computed in Quantum Electrodynamics 
(QED) - the most precise physical theory of the present time. QED is a highly accurate theory 
of electromagnetic particle interactions where predictions for physical observables are calculated 
within perturbation theory in the form of a series in the expansion parameter of the theory, the 
fine structure constant a with a numerical value of a^^ = 137.036 at zero momentum. The fact 
that a is rather small makes the expansions well convergent. The muon anomalous magnetic 
moment has been computed through the expansion in the small ratio of lepton masses and the 
expansion in the coupling constant up to the order a^. The most recent result reads [rUlllU lT^ 

a^(QED) = Ai + A2{m^/me) + A2{m^/m^) + A^{m^/me,m^/mr) (1) 

where for instance 

Ai = 0.5 (^^^ - 0.328478965... + 1.181241456... (^^y 

-1.5098(384) (^^y + 4.393(27) ■ 10^^^. 

The theoretical expressions for the coefficients of this series involve complicated multidimen- 
sional integrals the evaluation of which is the main difficulty in the calculation of a series 
expansion for the muon anomalous magnetic moment. A part of the coefficients in Eq. (P) 
have been calculated analytically whereas the numbers in the higher order terms have to be 
evaluated numerically [121 Ej . 

With the advent of powerful computers a number of problems can now be approached 
through direct numerical calculations. The trajectories of satellites, for instance, are computed 
through direct calculations with computers. There are also important applications in Quantum 
Field Theory (QFT) especially, in the theory of strong interaction formulated on the lattice. 
However, at present the achieved accuracy is not yet sufficient leaving plenty of room for 
analytical computations. In addition, calculations in perturbation theory are often formulated 
in such a way that a numerical approach is not efficient or is unreliable due to the accumulation 
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of rounding errors. The situation is quite analogous to the problem of long-time weather 
forecasts in meteorology where the direct numerical approaches are extremely unstable from the 
computational point of view. Rounding errors of the numerical evaluation of a large number of 
Feynman diagrams are such that the final accuracy is low because of huge cancellations between 
separate contributions of separate diagrams. Although an analytic evaluation is definitely 
preferable it cannot always be done due to technical difficulties in calculating the relevant 
Feynman integrals. 

In Quantum Field Theory only a very limited number of problems have a full exact solution. 
Perturbation theory is a general tool for investigating the realistic physical situations. This is 
especially true in the theory of strong interactions and in the Standard Model in general. The 
classical examples are the analysis of e~^e~ annihilation into hadrons and the analysis of tau de- 
cays where perturbation theory is used in its most advanced form for the calculation up to very 
high orders of the strong coupling constant ITT)] . In the latter case such an analysis has lead 
to the most accurate extraction of strong interaction parameters at low energies ^Tj |^ . 
The attempt to extrapolate the calculations to all orders for some limited subsets of diagrams 
or contributions are also quite useful [201 lUl 1221 12^1 • The coefficients of the expansions are 
given by complicated integrals which have a useful graphical representation in terms of Feyn- 
man diagrams. Graphical techniques are very efficient and convenient in the bookkeeping of 
the various loop contributions and are in wide use. In higher orders these graphs have many 
loops which are associated with the loop integrals. One speaks of multi-loop calculations for 
the evaluation of the coefficients of the perturbative expansion in QFT. Due to big advances in 
symbolic computing the most advanced symbolic programs can nowadays generate the terms of 
the perturbative expansion without actually drawing any graphs. Nevertheless, the terminology 
has survived. Nowadays the methods of QFT (and its terminology) is starting to penetrate 
into fields other than high energy physics such as polymer physics, hydro- and aerodynamics 
and weather forecasts which are much closer to real life. Therefore, the calculation of the com- 
plicated integrals associated with Feynman diagrams in a perturbative expansion is important 
in high energy physics as well as in many other areas of physics. 

While the problem of calculating multi-loop integrals will not be solved in a general form in 
the foreseeable future, there is a continuing interest in the analysis of particular diagrams in the 
efforts to enlarge the number of known analytical results (as a review, see e.g. [24]). Especially 
promising in this regard is the analysis of subclasses of diagrams with simpler topological 
structures but with an arbitrary number of loops. The most promising in this respect is the 
simple topology of the so-called generalized sunrise diagrams (sometimes also called sunset, 
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water melon, banana, or basketball diagrams). They will be referred to as the sunrise- type 
diagrams in this report. They can be computed for an arbitrary number of loops which makes 
them an ideal tool to obtain clues to the structure of more complicated multi-loop diagrams 

Recently there has been a renewal of interest in the calculation of diagrams of the sunrise- 
type topology for different numbers of loops involving different loop masses and/or external 
momenta of the diagram (see e.g. [2HI ESI EDI EI] ) • The classical (or genuine) two- loop sunrise 
diagram with different values of the internal masses has recently been studied in some detail 
using the traditional momentum space approach (see e.g. j32l OHl OH IHSl IHEj and references 
therein). In the present report we review the configuration space technique (also called x-space 
technique) for the calculation of sunrise-type diagrams. The configuration space technique is 
best suited for the topology of the diagram. It considerably reduces the complexity of the 
calculation and allows for new qualitative and quantitative results for this important particular 
class of Feynman integrals [23 EH EH EHl HOI lU] • Configuration space techniques can be used 
to verify known results obtained within other techniques both analytically and numerically jl2l 
HSl ISl 1^ and to investigate some general features of Feynman diagram calculation |^ IHU} 
EHl mi UHl HSl Eni EII- We list features of the x-space technique in turn with the aim to show 
how efficient the method is. 

The diagrams with sunrise-type topology form a subset of the general set of diagrams with 
a given number of loops. They appear in various specific physics applications: 

• multi-loop calculations in general EH] 

• the analysis of static properties of baryons using sum rule techniques in QCD [SH ESI EHl 

EZl EHl EHl EOl EH Ea EB] 

• the properties of glueballs extracted from the perturbative analysis of two-point correla- 
tors t64tl65] 

• the multi-particle phase space for multi-body decays of a single-particle phase space 

EH EH 

• lattice QCD calculations [HHl EO] 

• mixing of neutral mesons jTTj 

• Chiral perturbation theory (ChPT) and effective theories for Goldstone modes in higher 
orders of the momentum expansion [SZl EEl IZ2] 
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• analysis of exotic states in QFT: multi-quark states in QCD or pentaquarks using sum 
rule techniques |ZSllZlllZ3IZnilZZllZHllZnilSnilHI] 

• general questions of QFT: sum rules in two-dimensional QED |82L and properties of 
baryons in the large Nc limit of QCD jSH Ei] 

• effective potentials for symmetry breaking jHSl EH IHHl IHH IHOI 

• finite temperature calculations [MlEIlinSinSlinillMlEHlinZI 

• applications in nuclear physics [Ml E3 11001 IMl 11021 1103 

• phase space integrals for particles in jets where the momentum along the direction of the 
jet is fixed [TIHj 

• applications in solid state physics |lU5j 

The principal aim of this report is to assemble all the necessary tools needed for the com- 
putation of multi-loop Feynman diagrams with sunrise-type topology. We discuss the whole 
known spectrum of different methods to analyze these integrals. Amomg these are concise an- 
alytical evaluation techniques, expansion in small parameters such as masses and momenta (or 
inverse masses and momenta), expansions in special kinematical regions as the threshold regime 
in the Minkowskian domain, integral representations for integrals including their discontinuity 
across physical cuts in terms of analytic functions of the external momentum, and, finally, ef- 
ficient, fast, simple and stable numerical procedures. In the Euclidean domain the numerical 
procedures derived from our representation are efficient and reliable, i.e. stable against error 
accumulation. In order to acquaint the reader with the tools given given in the report we start 
with basic notations and relations. 

We present a comprehensive report on configuration space techniques including many tech- 
nical details such that an interested reader can use this report as a practical guide for practical 
calculations. We have included a great deal of mathematical material such that the report is 
self-contained to a large degree. In most cases it is not necessary to consult mathematical hand- 
books in order to understand the calculations. We also give a rather large sample of worked-out 
examples. The calculational methods used are rather well-suited for further development and 
can easily be tailored to the further specific needs of the potential user. Using the results of 
this paper one can create one's own software for an efficient evaluation of the many quantities 
of interest that can be extracted from the analysis of sunrise-type diagrams. 
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The report is organized as follows. In Sec. 2 we briefly summarize some general notions and 
their application to sunrise-type diagrams. In Sec. 2.1 we introduce the configuration space rep- 
resentation of sunrise-type diagrams and fix our notation. In Sec. 2.2 we specify what is meant 
by computing a sunrise-type diagram. In Sec. 2.3 we discuss the ultraviolet (UV) divergence 
structure of sunrise-type diagrams and present recipes to regularize the UV divergences either 
by subtraction or by dimensional regularization. In Sec. 2.4 we comment on the spectral density 
of sunrise-type diagrams and its connection to phase space. Sec. 3 is devoted to some explicit 
configuration space calculations involving both analytical and numerical methods the results 
of which are compared to previously known results where other calculational techniques have 
been used. One of the examples is the direct computation of the spectral density of sunrise-type 
diagrams without taking recourse to Fourier transforms. The results are important for multi- 
body phase space calculations. In Sec. 4 we discuss methods to find asymptotic expansions 
in different kinematical regimes of mass and/or momentum configurations in the Euclidean 
domain. Sec. 5 contains a generalization to non-standard propagators and non-scalar cases. 
In Sec. 6 we generalize the configuration space technique to other topologies. Sec. 7 contains 
our conclusions. Some of the lengthier formulas are relegated to the appendices, together with 
useful mathematical material about Bessel functions and Gegenbauer polynomials and a short 
treatise about cuts and discontinuities as they are occuring in the main text. 



2 Basic notions and relations for sunrise-type diagrams 

Sunrise-type diagrams are graphic representations of the n-loop two-point correlation functions 
in QFT with {n + 1) internal propagators connecting the initial and final vertex. The well- 
studied genuine sunrise diagram shown in Fig. ^a) is the leading order perturbative correction 
to the lowest order propagator in 0^-theory, i.e. it is a two-point two-loop diagram with three 
internal lines. This diagram emerges as a correction to the Higgs boson propagation in the 
Standard Model. It also naturally appears in some effective theory approaches to critical phe- 
nomena and studies of phase transitions in QFT. The corresponding leading order perturbative 
correction in 0^-theory is a one-loop diagram which can be considered as a oversimplified case 
of the prior example. 

The two-loop case is a standard starting point for the calculation of radiative corrections in 
QFT. It emerges in a huge variety of physical applications. Fortunately it can be analytically 
computed in any desired kinematical regime for any values of the relevant parameters. In this 
respect it does not present a real challenge as far as multi-loop calculations are concerned. We 
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Figure 1: Genuine sunrise (a) and general topology of the class of sunrise-type diagrams (b) 

shall, nevertheless, discuss the genuine sunrise diagram in some detail in order to illustrate 
the efficiency of our configuration space methods. We shall compare our results with known 
exact analytical results obtained with other techniques. This provides a mutual check on the 
correctness of the results. A straightforward generalization of this topology to the multi-loop 
case is a correction to the free propagator in 0"+^-theory that contains n loops and (n+1) 
internal lines (see Fig. ^b)). 

A subclass of the general sunrise-type diagram shown in Fig. ^b) is the case when the 
two external momenta vanish. This subclass is referred to as the subclass of vacuum bubble 
diagrams with sunrise-type topology which will be referred to as vacuum bubbles. They have 
a simpler structure than the true sunrise-type diagrams since the number of mass scales is 
reduced by one for the vacuum bubbles. We shall frequently return to the discussion of vacuum 
bubbles in the main text. 

2.1 Definitions and notation 

A n-loop sunrise-type diagram is a two-point correlation function with (n + 1) propagators 
connecting points x and y and as such is explicitly given by a product of simple propagators 
built from the basic two-point correlation function D{x,y,m), 

n+1 

U{x,y) = l[D{x,y,mi) (2) 

i=l 

(and/or their derivatives if necessary). The quantity D{x, y, m) may be the propagator of a free 
massive particle with mass m or a more general two-point correlation function. The sunrise- 
type diagrams are the leading order contribution to a two-point correlation function of two 
local currents jn+i (x) of the form 

jn+l (x) = 01 ■ ■ ■ X>^„+, (pn+l (3) 

where the fields (pi have masses and where is a derivative with multi-index = {/ii, . . . , fik} 
standing for = /dx^-^ . . -dx^^. The sunrise- type diagrams are contained in the leading 
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order expression for the polarization function 



^{x,y) = {Tjn+i{x)jn'+i{y)) (4) 

where the brackets mean quantum mechanical averaging over the ground state which is explicitly 
given by a product of propagators and/or their derivatives, 

U{x,y) =V^^^^{x,y,mi)- ■■Vf,„^^^^^^{x,y,mn+i). (5) 

As the standard case we will consider a translation invariant situation in which the propagator 
depends on the difference of the arguments only, D{x,y,m) = D{x — y,m). An exception to 
this important standard case would be the existence of an arbitrary external field. However, in 
the present report we will not further discuss this possibility. In the standard case one of the 
vertices can conveniently be placed at the origin, say y = 0. The propagator is then a function 
of X only, D{x,m). The basic expression for the sunrise-type diagram in configuration space 
reads 

n+l 

U{x) = l[Dix,mi). (6) 

i=l 

D{x,m) represents a free propagator of a massive particle with mass m in D-dimensional 
(Euclidean) space-time. It is given by 

where we have used D = 2X + 2. K\{z) is the modified Bessel function of the second kind 
(sometimes also known as McDonald function, see Appendix^ defined by Eq. ()A8|) . 

The massless propagator can be obtained from Eq. (jZj) by taking the limit m ^ (or more 
precisely the limit mx — at fixed x in terms of the dimensionless quantity mx). It reads 

and where r(A) is Euler's Gamma function. 

In some physics applications one may have more general basic two-point functions ( "modified 
propagators") for massive particles. An example is the calculation of the three-point function 
necessary for the determination of a particle form factor at zero momentum transfer when the 
external momentum of the current vanishes. This momentum configuration reduces the diagram 
to a sunrise- type diagram with more complicated basic two-point lines (see e.g. [106]). This 
phenomenon frequently occurs also in calculations of various quantities in Chiral Perturbation 
Theory (ChPT) |l()7j . Formally, these two-point functions cannot be considered as propagators 
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of a physical particle but emerge as effective basic two-point functions. The simplest modi- 
fication of the basic propagator which frequently appears in applications is the occurrence of 
higher powers of the standard propagator. In momentum space it has the form 

D^''\p,m) = ^^r-T- (9) 

For integer fi these cases appear if one e.g. considers mass derivatives of the propagator as they 
are needed in certain parameter expansions. For instance, one can consider 

^ (jp, _|_ ^2-) ^ ) 

In the same way a derivative of the momentum itself will also increase the power of the prop- 
agator. For instance, the second order scalar derivative □ = d^d^ of the standard propagator 
will result in 

^ + m? (p2 -)- m?Y {P^ + ^^)^ 
Higher powers of propagators occur in many applications. Therefore, any calculational tech- 
nique should be well suited to treat such higher powers. It turns out that the configuration 
space method very naturally accommodates higher order derivatives. 

In configuration space higher powers of the standard propagator with mass m are explicitly 
given by (cf. Eq. ((Zj)) 

1 f e*(P-^)rf^p 1 fm\^-^' 



the additional power of denominator factors (or, equivalently, the number of derivatives in 
masses) is denoted in the corresponding Feynman diagram by dots placed on the line (though 
even the power can be non-integer). Note that the modified propagator has the same functional 
form as the usual propagator. The main difference is the change of the value of the index of the 
Bessel function Kx{mx) Kx_^{mx). This fact unifies and simplifies the analysis of sunrise- 
type diagrams that contain more complicated basic propagators (analytical expressions that 
correspond to one line of the diagrams) as they emerge in some physical applications or are 
used for the calculation of other diagrams through recurrence relations techniques 



2.2 Momentum versus configuration space representation: 

Fourier transform in the evaluation of sunrise— type integrals 

The calculation of loop diagrams in perturbative Quantum Field Theory is associated with 
the evaluation of complicated multi-dimensional integrals, see e.g. jl()9j . It is obvious that the 
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expression for the sunrise-type diagram in configuration space given by Eq. © contains no 
integration at all even if it represents a multi-loop diagram. In this respect it is a kind of 
tree-level diagram that usually appears in the lowest orders of perturbative QFT. However, 
"to calculate" a sunrise diagram means in most cases to find certain integrals of Eq. over 
X with some weight functions depending on the quantity in question. The simplest situation 
corresponds to a computation of the values of the vacuum bubbles, i.e. expressions for the 
diagrams without external momenta (or with vanishing external momenta, see e.g. |ll()j ). In 
case of the sunrise-type topology the value of the vacuum bubble is given by integrals with a 
simple weight function. 



This integral provides the simplest example of calculating a sunrise-type diagram. We like 
to emphasize that this calculation is almost always a must as it allows one to determine the 
counterterms related to UV divergences of sunrise-type diagrams that always occur in realistic 
physical models. Due to rotation invariance the integral is one-dimensional as the angular 
integration can be explicitly done. Thus 



The possibility to compute all vacuum bubbles in terms of one-dimensional integrals for any 
number of loops is a big advantage of the configuration space technique. This guarantees 
that the complete renormalization program can be explicitly accomplished in terms of one- 
dimensional integrals. It should be clear that this is a huge simplification. 

The most familiar example for the evaluation of a sunrise-type diagram is the procedure of 
finding its value in momentum space, or put differently, the evaluation of its Fourier transform. 
It corresponds to non-bubble diagrams with two external lines as the integrals then requires 
the computation of the Fourier transform of n(x). 



The required integrals are basically scalar integrals. This makes the angular integration in 
Eq. ()13|) simple in D-dimensional space-time. The needed formula is the integration over the 
angular variables d^x which, in explicit form, reads 




(13) 




(14) 




(15) 




(16) 
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where p = \p\, x = \x\, and J\{z) is the Bessel function of the first kind. The integration 
/ d^x over the rotationally invariant measure d^x on the unit sphere in D-dimensional (Eu- 
chdean) space-time means integration over angles only. The generalization of Eq. (fT^ to more 
complicated integrands with an additional tensor structure x^^ ■ ■ ■ x^'' is straightforward and 
leads to different orders Jx^i{z) of the Bessel function for different irreducible tensors (denoted 
by the number /) after angular averaging (integration over angles). The corresponding order 
of the Bessel function can be inferred from the expansion of the plane wave function e**-^'^-' 
which appears as a weight function upon computing the Fourier transform in a series of Gegen- 
bauer polynomials C^{w) (see Appendix |B}. The Gegenbauer polynomials are orthogonal on 
a D-dimensional unit sphere, and the expansion of the plane wave in terms of Gegenbauer 
polynomials reads 

^.(p..) ^ r(^) ^P^y'±^^iX + l)J,^,ipx)C^ (^^] . (17) 

Besides the definition of the Gegenbauer polynomials by a characteristic polynomial as given 
in Eq. (jBljl . the expansion in Eq. ()17|1 can also serve as a definition of these polynomials. 

The representation in Eq. ()17|) allows one to single out an irreducible tensorial structure 
from the angular integration given by the Fourier integral in Eq. (llHj) . Integration techniques 
involving Gegenbauer polynomials for the computation of multi-loop massless diagrams are 
described in detail in jlllj where many useful relations can be found (see also |112j where the 
calculation in momentum space was considered). One arrives at a representation of the Fourier 
transform of a sunrise- type diagram that is given in terms of the one-dimensional integral 

Il{p) = 2vr^^'^"(Y) ^^A(pa;)n(x)x'^+Mx 

= ^^^^^ [y) ^'^A(pa;)D(x,mi)---D(x,m„+i)x2^+Mx (18) 
which is simple to analyze. The representation given by Eq. (jTH)) is quite universal. 

2.3 Regularization and renormalization 

While the configuration space expression for the sunrise-type diagram contains no integration 
it does not always represent an integrable function of x in configuration space for general 
values of the space-time dimension D. The expression given by Eq. (jHI) can have non-integrable 
singularities at small x for a sufficiently large number of propagators if the space-time dimension 
is D > 2. 
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Upon multiplication of many propagators one may get a function which is singular at small 
X. Each propagator by itself can be considered as a distribution. However, distributions do not 
form an algebra and their multiplication is not well-defined jll3l . In this respect the polarization 
function 11 (x) in Eq. (jHl) is not completely defined as a proper distribution. In attempting to 
integrate such a function over the whole x-space we encounter infinities in the form of ultraviolet 
(UV) divergences |114j . Therefore, the computation of the Fourier transform of n(x) requires 
regularization (for instance, dimensional regularization) and subtraction. Configuration space 
provides a nice environment for this. 

Note that all these statements depend on the dimension of space-time. The vicinity of 
D = 2 is special since the leading singularity of the basic propagator is only logarithmic and the 
multiplication of propagators does not lead to non-integrable singularities. Therefore, the UV 
properties of two-dimensional models that are often considered as simplified models in QFT or 
that emerge in real physical applications as a result of some given approximations |115|lll6|lll7j . 
They are much softer than those of realistic four- dimensional models (see e.g. |118j ). 

In general, however, one has to consider an arbitrary space-time dimension D. The modern 
way of dealing with singularities in QFT in particular in multi-loop calculations is mostly based 
on dimensional regularization |119|ll20]ll21j (as a review see |122j ). Dimensional regularization 
is characterized by considering the space-time dimension D to be an arbitrary parameter, for 
instance D = Dq — 2e where Dq is an integer dimension.^ In case of non-integer space-time 
dimension D the integration is well-defined, while the limit D ^ Dq is singular, consisting 
usually of poles in D — Dq. These pole parts are extremely simple to handle in configuration 
space. Indeed, in configuration space the regularization is given by adding the (5-function 6{x) 
and its derivatives to n(a;). As we will see in the following, this procedure corresponds to adding 
a polynomial in momentum space (local counterterms according to the Bogoliubov-Parasiuk- 
Hepp-Zimmermann (BPHZ) theorem j ll3t I123t I124t I125j ). The finite parts can be found 
numerically for any number of loops by doing a single one- dimensional integration. Therefore, 
the renormalization in x-space becomes 

n^(x) = U{x) + CoS{x) + Cin6{x) + . . . + Crn''6{x) (19) 

where □ = d^d'^. The coefficients Ci are functions of the regularization parameter e which are 

singular in the limit £ — > 0. The extraction of the coefficients Ci is explained in the subsequent 

^The usual choice of the integer space-time dimension Dq for real physical applications is Dq ~ 4, but other 
space-time dimensions Dq — 2, Dq ~ 3, Dq = 5 or Dq = 6 are being used in other applications. Dq = 3 is 
natural for near threshold expansions or expansions in the Minkowskian domain in general. 
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sections for some special cases. The number r of counterterms is determined by dimensional 
arguments such as power counting and can be easily related to the number of propagators and 
the space-time dimension. The actual number of poles depends on the mass configuration. 



2.4 The spectrum as discontinuity across the physical cut 

From the physics point of view one important aspect of our analysis of sunrise-type diagrams is 
the construction of the spectral decomposition of the diagrams. For the two-point correlation 
function we determine the discontinuity across the physical cut in the complex plane of the 
squared momentum, = —rn? ± iO which is referred to as the spectral density 

p{s) = -^Discn(|}) 

ZTTZ p^=-s 

(cf. Appendix in)). Note that the spectral density of a sunrise-type diagram is finite for any 
number of loops. It turns out that the configuration space technique allows one to compute 
the spectral density in a very efficient manner. The analytic structure of the correlator n(x) 
(or the spectral density of the corresponding polarization operator) can be determined directly 
in configuration space without having to compute its Fourier transform first. The technique 
for the direct construction of the spectral density of sunrise-type diagrams introduced in [37] 
is based on an integral transform in configuration space which in turn is given by the inversion 
of the relevant dispersion relation. 

The dispersion representation (or spectral decomposition) of the polarization function in 
configuration space has the form 

n(x) = / p{s)D{x,y/s)ds = p{m^)D{x,m)dm^ (20) 

JO JO 

where ^/s = m. This representation was used for sum rule applications in jl2(i| I127j where 
the spectral density for the two-loop sunrise diagram was given in two-dimensional space- 
time |128j . The representation in momentum space is more familiar and is referred to as the 
Kallen-Lehmann representation jl29j of the two-point correlation function. In the Euclidean 
domain it is given by 

JO s 

This expression can of course be obtained by taking the Fourier transform of both sides of 
Eq. (1201). 

With the explicit form of the propagator in configuration space given by Eq. (0), the repre- 
sentation in Eq. ()20p turns out to be a particular example of the Hankel transform, namely the 
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/^-transform |13U1 1131j . Up to inessential factors of x and m, Eq. (j^UI) reduces to the generic 
form of the i^-transform for a conjugate pair of functions / and g, 

9{x) = / f{y)K^{xy)y/xydy. (21) 
Jo 

The inverse of this transform is known to be given by 

]^ rc+ioo 

fiy) = — g{x)I^{xy)y/xydx (22) 

TTZ J c—ioo 

where Iu{x) is a modified Bessel function of the first kind and the integration runs along a 
vertical contour in the complex plane to the right of the right-most singularity of the function 
g{x) |131j . In order to obtain a representation for the spectral density p{rn?) of a sunrise-type 
diagram in general D-dimensional space-time one needs to apply the inverse i^-transform to 
the particular case given by Eq. (j20|) . One has 

m^p{m^) = ^—^ / Il{x)x^+^Iximx)dx. (23) 

I J c—ioo 

From Eq. ()23p we obtain an explicit analytical representation for the spectral density p{s) as a 
contour integral of the polarization function which reads 

/'(^) = / h{xV^Mx)x^+^dx (24) 

IS ' J c—ioo 

where Ix{z) is a modified Bessel function of the first kind and the integration runs along a 
vertical contour in the complex plane to the right of the right-most singularity of n(x). 

The inverse transform given by Eq. completely solves the problem of determining 

the spectral density p(s) of the general class of sunrise-type diagrams by reducing it to the 
computation of a one-dimensional integral along the contour in a complex plane. This is valid 
for the general class of sunrise-type diagrams with any number of internal lines and different 
masses. For {n + 1) internal lines with (n + 1) equal masses m the spectral density reads 

™A(ra+l) pc+ioo 

P(^) = ,(2^)nA...i,V2 l_ ^A(a: V^) iKximx)r' x'^^^dx. (25) 

Because the contour can bypass the area of small values of x, the integral is finite, as it was stated 
before. The spectral density in turn can be used to restore the finite part of the correlation 
function in the momentum space representation. The path traced here is is an alternative to 
the calculation of the Fourier transform for sunrise-type diagrams. 

In the standard, or momentum representation, the polarization function Il{p) is calculated 
from a n-loop diagram with n Z)-dimensional integrations over the entangled loop momenta. 
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Figure 2: The spectral density for the four-hne sunrise-type diagram with equal masses for 
D = e = 2.718 . . ., D = 3, and = vr = 3.14 . . . space-time dimensions. 



It is clear that the momentum space evaluation becomes very difficult when the number of 
internal lines becomes large. 

In order to demonstrate the applicability of our method described in more detail in [SZl EHl 
EH], we show in Fig. |21the result for the spectral density for the four-line sunrise- type diagram 
with D = e = 2.718 . . ., D = 3, and D = tt = 3.14 . . .. We have chosen these exotic values for D 
to demonstrate the power of configuration space techniques. To the best of our knowledge there 
is no other method which allows one to compute the spectral density with such completeness 
and ease, i.e. through a one-dimensional integral. The technique was first suggested in |37, and 
has been studied in detail in Refs. [HHl 1^. 

Note, finally, that the spectral density of a sunrise-type diagram is a representation of 
multi-particle phase space with the number of particles equal to the number of lines of the 
diagram. Our formula for the spectral density is a striking example of how the configuration- 
space technique can be used to determine the phase space of physical processes. Efficient 
techniques for the calculation of the multi-particle phase space are important for Monte-Carlo 
simulations and cross section calculations for multi-particle final states as they will occur at 
the upcoming Large Hadron Collider (LHC) under construction at CERN. 
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3 "Tool box" for a practical analysis 



In this section we give examples and practical prescriptions of how to proceed with the calcu- 
lation of sunrise-type diagrams containing only a single Bessel function and powers of x. The 
massless case can be solved analytically in closed form for any number (ra + 1) of internal lines 
since it contains only a single Bessel function and powers of x. The case for one massive line 
and n massless lines (with two Bessel functions) is completely solvable analytically as well. 
The first nontrivial situation emerges with two massive lines. In this case a complete analytical 
solution is unknown in most of the cases, the exception being the case of equal masses that 
can be integrated analytical through hypergeometric functions. However, when one separates 
the diagram into a singular or pole part and a finite part, both parts can be calculated either 
analytically (for the pole part) or numerically (for the finite part) in very general situations. 
This will be demonstrated in this section. 



3.1 Extraction of poles in dimensional regular izat ion 

Renormalization counterterms can be constructed analytically for general mass configurations 
and any number of loops. As mentioned earlier, in dimensional regularization the UV diver- 
gences become manifest as poles in e ^ D — Dq. In configuration space the UV divergences 
are related to short distances. In order to analyze the structure of singularities at small x, one 
has to expand the massive propagators at small x (or, effectively, in the masses). In order to 
obtain the requisite conterterms one then integrates over x which is analytically feasible for 
D 7^ Dq, i.e. e 0. Nevertheless, one should ensure convergence also at large x which requires 
an appropriate IR regularization. Only this guarantees that singularities that emerge after 
integration are related to small x divergences only (see e.g. |132j ). This technical constraint 
can be met for instance by keeping one (or two) of the massive propagators in the integrand 
unexpanded. Since the massive propagator falls off exponentially at large x, this procedure is 
sufficient to ensure convergence at large x. 

The possibility of retaining one propagator unexpanded exists because the integrals of one 
or a product of two Bessel functions are known analytically |1331 11341 1135j . This nice feature 
provides a tool for obtaining counterterms analytically. The necessary formulae read 



x^-^K,{mx)dx = 2'^-2m-'^r (^^^) r (^^i^) (26) 

22a-3 

m2"r(2a) 



/•oo 22"-3 

/ x^"-^K^{mx)K^{mx)dx = -———T {a + fi)r {a)r {a)r {a - fi) . (27) 
Jo m'^°'i 2a 
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In general, the finite parts of tlie relevant diagrams resulting from the subtraction of coun- 
terterms can be calculated only numerically. Still there are some simple examples where an 
analytic calculation is feasible. One of these examples is the three-loop sunrise-type diagram 
with two massive and two massless lines at vanishing external momentum, the three-loop vac- 
uum bubble. There exist an analytical expression for the value of this particular vacuum bubble 
in configuration space. It is given by 

n(0) mm^, Of,^. = / ( "":':fj''T> l ^ ( ^] ' (28, 



While the angular integration in D-dimensional space-time is trivial, the problem of the residual 
radial integration is solved by using Eq. ()27|) . The result for the integral in Eq. (j^H|) reads 

f,,r,. ^ irn'r^' r(A)^r(l-A)r(l-2A)^r(l-3A) 
^ ' (47r)3{A+i) r(A + i)r(2-4A) ■ ^ ^ 



This result corresponds to the quantity Mi in Ref. |136j where it constitutes the simplest basis 
element for the computation of massive three-loop diagrams in a general three-loop topology. 
Again, any number of massless lines can be added. 

The result given by Eq. ()29|) is obtained in a concise form and valid for any dimension D. 
Therefore its pole parts can be explicitly singled out through a direct expansion into a Laurent 
series near an integer value of the space-time dimension D. Taking D = A — 2e (i.e. A = 1 — £:) 
and expanding the result for small e one readily finds 

- (m^-^- ( 1 7 25 + 6C(2) 5 - 42C(2) - 56C(3) \ 

^ ' 3(47r)3(2--) + 2£2 + 4£ 8 ^ ^ ') ' 

Using this example one can take a closer look at the two-dimensional case. As was stated 
before, the case Dq = 2 is much softer and the UV divergences are practically absent. Hence 
one could come to the conclusion that the expression given by Eq. should be perfectly finite 
a.t D = Dq = 2. In terms of D = 2A + 2 the limit corresponds to taking A ^ 0. Therefore, one 
might be surprised to find a singularity r(A)^ in Eq. ()29|) . However, this is not a UV singularity. 
Because massless propagators are not defined in the two-dimensional world, the limit A — is 
ill-defined in this case. Indeed, looking at Eq. (jH)), the singularity of this sunrise- type diagram 
is exactly the product of two singular massless propagators. Power counting for the integral 
expression in Eq. (jH)) shows that the singularity is caused by the infrared (IR) region. As this is 
an IR problem, the afore mentioned BPHZ theorem does not apply and the divergence cannot 
be treated with local counterterms in the usual manner. 

Returning to the more general case, a closed- form evaluation (as in Eq. (|29|) ) is not possible 
and an explicit subtraction is required to separate singular and finite parts. In our case it is 
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convenient to use momentum subtraction which in fact is the oldest renormahzation method. 
Momentum subtraction consists of subtracting a polynomial at some fixed momentum point. 
For massive diagrams the expansion point p = is safe (i.e. has no infrared (IR) singularity), 
and the prescription is realized by expanding the function 

[^y'jxipx] (30) 

(which is the kernel or weight function of the integral transformation in Eq. ()18|)) in a Taylor 
series around p = in terms of a polynomial series in (cf. Eq. (IA4jl in Appendix El- The 
subtraction at order is achieved by writing 

and by keeping A^ terms in the expansion on the right hand side. After substituting the 
expansion Eq. ()31|) into Eq. (fTHj) one obtains the momentum subtracted polarization function 

n..om(p) = %)-E^(^j n(p)lp^=o (32) 

which is finite if the number of subtractions A^ is sufficiently high. The function Il{p) is divergent 
as well as any of the derivatives on the right hand side of Eq. (j32p . The divergences require 
regularization. However, the difference, i.e. the quantity Hmomip) is finite and independent of 
any regularization used to give a meaning to each individual term in Eq. ()32|1 . Note that the 
expansion in Eq. (jHTj) is a polynomial in in accordance with the general structure of the 
i?-operation |114j . The number A^ of necessary subtractions is determined by the divergence 
index of the diagram and can be found according to the standard power counting rules. The 
subtraction at the origin p = is allowed if there is at least one massive line in the diagram 
along with an arbitrary number of massless lines. If there are no massive internal lines at all, 
the corresponding diagram can easily be calculated analytically and the problem of subtraction 
is trivial. After having performed the requisite subtraction, one can take the limit D Dq 
in Eq. (fTHj) where Dq is an integer. The diagram as a whole becomes finite after subtraction 
which refiects the topology of the sunrise-type diagram: there is no divergent subdiagram in 
the sense of Bogoliubov. 

In selecting the momentum subtraction in order to regularize the singularities we have 
selected a particular subtraction or renormahzation scheme. The coefficients Ci in Eq. (fT^ 
depend on this choice. While the highest coefficient related to the strongest singularity is 
unique, the other coefficients are scheme dependent. The standard scheme used in the literature 
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^) J\{px)-Y. 



t^, mx + k+i) 
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Y 
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(31) 
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is the MS scheme. The transition from the momentum subtraction scheme to the MS scheme 
is achieved by the famihar redefinition of the measure 

d^x {Ann^e-^^fd^x (33) 

with /i being the renormahzation scale within dimensional regularization and 7^; being Euler's 
constant •je = 0.577 . . . 



3.2 Singular and finite parts for particular parameter values of the 
genuine sunrise: warm-up example demonstrating the technique 

As an example of calculating the counterterms and evaluating the remaining finite part of the 
diagram we present the special case of the genuine sunrise diagram with three massive lines, 
i.e. the two-loop sunrise-type diagram. For special values of the external momenta a closed- 
form expression for the genuine sunrise diagram is known within dimensional regularization 
using momentum space techniques [^3]. We shall reproduce this result in the following. The 
one-dimensional integral to be analyzed is given by 



00 / pjr 
V 2 



n(p) = 2tt^+^ I { — ) Jxipx)D{x, mi)D{x, m2)D{x, m3)x^^+^dx. 



(34) 



To determine the finite part we first use momentum subtraction and split Il{p) up into its finite 
and infinite (but dimensionally regularized) part, 



U{p) = n 

mom 

(p) + nsing(p) 



(35) 



where tlmomip) is a momentum subtracted polarization function and Ilsing(p) is a counterterm 
in dimensional regularization. Power counting shows that only two subtractions are necessary. 
The explicit expression for the momentum subtracted polarization function reads 



2n 



A+1 



Y ) -Mpx) 



D{x, ■mi)D{x, m2)-D(x, ■m3)x^^'^^dx 



J 1 



27r 



A+1 



px 

Y 



Jx{px) 



r(A + i) 4 r(A + 2) 



xD(x, mi)D{x, m2)-D(x, m'i)x^^^^ dx 



(36) 



while the pole part is given by a first order polynomial in p^ of the form A + Bp^ . The actual 
values of the coefficients A and B depend on the regularization scheme that has been used. We 



23 



use dimensional regularization and obtain 



nsing(p) 



A + p^B 



2^A+1 ^oo 




4r(A + 2) Jo 



oo 



x'^D{x, 'mi)D{x, m2)-D(x, ms)^' 



(37) 



With the help of the representation given in Eq. (jHTj) we can apply our strategy in a straight- 
forward manner. The coefficients A and B are simple numbers independent of p. They contain 
divergent parts (regularized within dimensional regularization) and need to be computed only 
once to recover the function Il{p) for any p. In the momentum subtracted part one can forego 
the regularization (it is finite by the /^-operation) and perform the one- dimensional integration 
numerically ioT D = 4 (i.e. A = 1). 

In the particular case of the genuine sunrise diagram the necessary integrals are known 
analytically in closed form and can be found in integral tables (see e.g. |134j ). However, using 
these tables may not always be convenient and we therefore again present a simplified approach 
which allows one to deal with even more complicated cases in a simpler manner. Let us specify 
to the particular case p = mi + m2 — rris (which corresponds to a pseudothreshold) where an 
analytical answer for the integral exists [33]. For simplicity we choose mi = m2 = m,3/2 = m. 
Then p = and nmom(O) = (which is a regular function at this Euclidean point). In this 
case the counterterm Bp^ vanishes because the quantity B is finite at finite e. We thus only 
have to consider the coefficient A. On the other hand our considerations are completely general 
since for arbitrary p one only requires more terms in the p^-expansion. For the special mass 
configuration considered here the result of [3^ reads 



This result can be extracted from the first line of Eq. ()37p using the integral tables given in |134j . 
However, in the general mass case the necessary formulas are rather cumbersome. Even for the 
special mass configuration considered here they are not so simple. We therefore discuss a short 
cut which allows one to obtain results immediately without having to resort to integral tables. 
What we really need is an expansion in e. Basically we need the integral 



m"^(0) = n 



A-2e 



m2-4^r2(l +e) r 3 8 In 2 
{l-e){l-2e) ^ 






(39) 



which is of the general form 




(40) 
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{n = X and p = 1 — A in our case). For a product of two Bessel functions in the integrand 
(without the last one in the above equation) the resuh of the integration is known analytically 
and is given by Eq. ()27|). Let us reduce the problem at hand to Eq. ()27|) and do our numerical 
evaluations with functions in four-dimensional space-time where no regularization is necessary. 
To do so we subtract the leading singularities at small ^ from the last Bessel function in 
Eq. (j^0|) using the series expansion near the origin given by Eq. (jA10|l . After decomposing the 
whole answer into finite and singular parts according to (27r)^'^y4 = F + S (where the total 
normalization of has been adopted in the definition of F and S) we find for the singular 
part 



S 



r(A + i) 

(27r)^m2-4- 

r(A + i) 



x-^i^-^)-^Kx{mx)Kx{mx) 



r(A) 



1 + 

e 



[mx) 
1- A 

rfi 



i-A r(i + A) 



r(i-A) 
r(i + A) 



(mx) 



2A 
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e 



Oie). 



(41) 



{l-e){l-2e) 

The pole contributions coincide with the result in Eq. ()38|1 while the finite part is different. It 
is corrected by the finite expression 

(27r)^m2^ 



r(A + l) 

X < {nixYKx{2mx) 



x^(^-^)-^Kx{mx)Kx{mx) 

r(A) - 



[mx] 



rri 



i-A r(i + A) 



dx. 



Because F is finite (no strong singularity at small x) one can put A = 1 to obtain 

1 r 



OO fjqf 

—Ki{x)Ki{x) {xKi{2x) 

X 



l + x^(-l + 27B + 21nx) 



Doing the numerical integration in Eq. 



F = levr^m^ 



results in the expression 
0.306853... 1. 



(42) 



(43) 



(44) 



One can restore all e-dependence in the normalization factors as in Eqs. and ()4H1 because 
F is not singular in e, so that this change of normalization is absorbed by the 0{e) symbol. 
One obtains 



TT 



4-2e 



m 



2-4£f2 



r^(i + ^) 



-0.306853, 



+ 0{e) 



(45) 



(l_e)(l-2£) 

Adding F from Eq. ()45|) to S from Eq. (PT|) one obtains a numerical result for the finite part 
of Eq. In order to establish the full coincidence with the result of Eq. ((211), ^ote that 



-0.306853 . . . 



;i -ln2) X 1.0000.... 



25 



Of course we have used the analytical expression (1— In 2) for the numerical quantity 0.306853 . . . 
for illustrative reasons because we knew the final answer. 

We emphasize that there is nothing new in computing the polarization function related to 
this diagram at any p. Some finite part appears from the momentum subtracted polarization 
function. One then needs another counterterm for non-zero p^. The computation of these 
additional terms proceeds in analogy to what has been done before. It is even simpler because 
the singularity at small x is weaker and only one subtraction of the Bessel function is necessary. 
This computational technique works for any complex value p^. 



3.3 Three-loop sunrise-type diagram for special momenta 

In a similar manner one can treat the three-loop sunrise-type diagram with two equal masses 
M, one mass m with values between zero and M and a massless line for outer momentum 
= — m^. This particular sunrise-type diagram was analyzed by Mastrolia [23 whose results 
we confirm. The starting integral expression reads 

n(p) = '^^^^^ [y] ^^A(pa;)^(2;,M)D(x,M)D(a;,m)D(a;,0)a;2^+^(ix 
2(27r)3(A+i) io 

The integral has the same form as in the previous example. It is obvious that we can proceed in 
the same way as before: first expand the integrand in p and m in order to find the counterterms, 
and second determine the finite contribution for p"^ = —m^ numerically. However, in this case 
we have to expand the Bessel function J\[px) up to the second order in p^. For the singular 
part we obtain 

271-^+1 roo 

— r / D[x,mi)D[x,m2)D[x,m^)x 
1) 



(^y) ^ Mpx)Kx{Mx)Kx{Mx)Kx{mx)x^~^^dx. (46) 
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dx 
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r(A + 

27]-^+l roo 

-p^ I o\ / x^D{x,mi)D{x,m2)D{x,m'i)x'^^^^dx 

41 (A + I) Jo 

, n\ / x'^D{x,mi)D{x,m2)D{x,m3)x'^^'^^dx (47) 
321 (A + 3j JO 



(47r)3(2-^)r(2 - 2£)r(2 - e) I 6^3 + 24^2 
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+\ (-^(2 + y^) + ^(1 + 2|/2)C(2) - y2 In y + In ^ + ^2 ^ 
^ _ 2M2-6^ /2 + I/2 ^ 2 + ^/2_6'^/2lny^^ 



(47r)3(2-£)r(2 - 2£)r(3 -e) \ lle^ lie 

2M-^' / 1 ^ 

(47r)3(2-^)r(2 - 2e)T{A - e) [Ve ^ 



/ 1 \ 



and 2/ = m/M. 



3.4 Infinite parts of sunrise- type diagrams with different masses at 
any number of loops: general techniques 

If one deals with sunrise-type diagrams with different masses one has to be careful to use a 
symmetric subtraction procedure. Otherwise one looses the symmetry of the original diagram. 
Such an asymmetry would arise if one would expand one of the Bessel functions in the integrand 
of the sunrise-type diagram as was done in the previous example for the massive line with mass 
m. If one wants to obtain a result which is symmetric under the exchange of the masses m, 
as expected from the topology of the diagram and from the initial form of the integral we 
have to subtract a counterterm which is symmetric in these masses. In calculating the genuine 

2 2 

sunrise diagram one can follow this procedure by using a damping function such as ^ to 
suppress contributions for large values of x, corresponding to IR singularities. This is done by 

2 2 2 2 

introducing the factor 1 = ^ ^ into the integrand. In this case all massive propagators 

2 2 

along with the factor ^ can be expanded in terms of small x values. The expansion of the 
propagators is straightforward and can be obtained from Eqs. ()A4|1 and (jAlOjl in Appendix 1X1 
The final integration can then be done formally by using the identity 

roo „ 1 

/ x^-^e"^ ^ dx = -fx-PT{p/2) (49) 
Jo 2 

which is nothing but the definition of Euler's Gamma function. 

As an illustration we compute the four-loop bubble diagram which was also calculated by 
Laporta |137j whose results we confirm. In the sample results below we extract the normaliza- 
tion factor 

for the n-loop (or {n + l)-line) sunrise- type diagram. Extraction of this factor renders all pole 
parts to be rational numbers. The results of the calculation are 

1 



AfiU,{p') = - + 0{e 
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MMp' = 0) = „.{ 5 _ 35 4565 _ 583451 

^ 2 2^ ef 5 45 4255 1061471 ^, q, 

/ 2 sfS 33 1247 180967 8985171 

/ 2 2n 8f 3 16 49 6967 17060631 



Ar6n6(p' = 0) 



_ 7 _ 133 _ 238 77329 
"2^6 " 6^5 ^ 3^1 ~ 360^3 
21221921 2596372387 
43200^2 2592000e 



0{e^'). (51) 



The coefficient of the leading singularity in e is independent of p^. In Appendix |Dl we list 
results for unequal mass configurations up to four-loop order. When setting all masses equal 
the results of the general mass case can be seen to agree with the above equal mass results. 

3.5 Finite part for the genuine sunrise with three different masses 

In configuration space the genuine sunrise diagram with three different masses mi, m2, and 
is given by 

Il{x) = D{x,mi)D{x,m2)D{x,m^). (52) 
The Fourier transform of this polarization function in the Euclidean domain reads 

II(p) = 27r^+^ (^y ^ ^ Jx{px)D{x, mi)D{x, m2)D{x, ms)x^^+^dx. (53) 
In dimensional regularization the singular parts [D = 4 — 2e, i.e. X = 1 — e) are given by 

fr ^ N _ /^"^' [ ml + ml + ml ^ 1 mf ln(m^e^V2/i) \ \ 

lining W - 1 5^2£2 + s 128 1024 J j ' 

where 7' = 7£;/2 — 3/4. The counterterms for the Fourier transform are obtained as explained 

2 2 2 2 

before by introducing 1 = e'^ ^ ^ and expanding the relevant Bessel functions together with 

2 2 

X _ xhis is the reason for the explicit dependence on the arbitrary parameter /i. In this way 
one readily finds the pole parts for any diagram with any mass configuration, i.e. one finds the 
coefficients Ci of the counterterms in Eq. ()19j) . Even if an external momentum goes through 
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the diagram, the pole parts can be still calculated easily since they are given by integrals over 
the expanded integrand. The momentum subtracted finite part is given by 

D{x, mi)D{x, m2)D{x, majx^^+^e'^'^'e-^'^'cix (55) 

N 

where is the order of the expansion needed to ensure integrability. Naturally has the same 
order as the one used earlier in order to extract the counterterms. It is obvious, though, that 
the expansion can terminate at positive powers of x. Finally, one can set e = to calculate the 
finite part numerically. As a technical remark note that most numerical integration routines 
will run into problems if the upper limit extends to infinity. However, already for values of 
the order x ~ 1 the integrand is negligibly small, so that the integration can be cut off at 
this point. Also the region around the origin might cause trouble for the convergence of the 
numerical integration. In this case the integration interval can be subdivided with an interval 
close to the origin, and the whole integrand in this subinterval can be expanded in x and 
integrated analytically. For specific values of p^, mi, 1712, and we were able to reproduce 
results given in the literature (see e.g. results given in Refs. I138j ). 

3.6 Calculation of the finite part: further examples 

The examples presented so far are well-known and have been obtained before by using tech- 
niques different from the ones presented here. While we can numerically compute any sunrise- 
type diagram with any arbitrary number of internal massive lines it is not always possible to 
find corresponding analytical expressions to compare with. Beyond two loops only a few exam- 
ples can be found in the literature. For instance, we have calculated the result for the vacuum 
bubble diagram shown in Fig. El (see also Sec. 15. 4p . For p we obtain (cf. Eq. ()A4|) ) 

(f)"''^'^^'^f(ATT)- <^^' 

One therefore has 

n(0) = — -/ (D^^\x,m)) D{x,m)x'^^+^dx (57) 

r(A -I- 1) JO ^ ^ 

where D^^\x,m) is defined in Eq. ()12j) . In order to compare with the literature we 
multiply with a relative factor A/sm^"^^^ (cf. Eq. fl50|) ) to obtain 

- Ar3m--II(0) ^ .p.^ .3^(1-.) /o ^^.(-)^-^(-)-^^^^^- (58) 



nfin(p) = 27r 



A+1 
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Figure 3: Three-loop massive sunrise-type bubble. A dot on a line corresponds to one differen- 
tiation with respect to the propagator mass squared 

The integral remains finite even for e — >■ and, after finite integration, we obtain the 
numerical result 

POD 

Bo = 4 Kl{x)Ki{x)x^dx = 2.1035995805 . . . (59) 
Jo 

A comparison with the results presented in Ref. |13(ij shows full numerical agreement. Since we 
know that the quantity Bq contains only a single transcendental number, namely C(3), where 
({z) is Riemann's ^-function jl3(i| 1139] . we can take its numerical value ({3) = 1.202056903 . . . 
to obtain the result 

5o = C(3)- 1.7500000... = ^C(3) 

which is nothing but the analytical result given in Ref. |136j . 

We conclude that the configuration space method is simple and allows one to numerically 
compute diagrams with sunrise-type topology in a very efficient manner. When the structure of 
the transcendentality of the result is known for a particular diagram (as in the latest example 
where only ^(3) is present) or when it can be obtained from an educated guess looking at 
the topology of the diagram |140j our numerical technique can be used to restore the rational 
coefficients of these transcendentalities which can be thought of as elements of the basis for a 
class of diagrams. 

3.7 Higher order expansion in e: 

general techniques and the four-loop vacuum bubble 

After having determined the singular parts of the Laurent series expansion using the damping 
factor method what remains to be done is to calculate the coefficients of the positive powers 
of e in the e-expansion including the finite term. Higher order terms in the e-expansion are 
needed if the sunrise-type diagram is inserted into a divergent diagram or when one is using 
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the integration-by-parts recurrence relation which can generate inverse powers of e. In order to 
determine these higher order terms one needs to resort to the Bessel function method. What 
is technically needed is to develop a procedure for the e-expansion of Bessel functions. 

The 0{e) term can be obtained analytically. Indeed, within dimensional regularization the 
propagator in configuration space contains Bessel functions with non-integer index depending 
on the regularization parameter e, as they occurred for instance in Sec. 13.61 at an intermediate 
step. In order to expand the Bessel function in the parameter e entering its index, we use 
Eq. ()A13|) for the derivative of the modified Bessel function of the second kind with respect to 
its index near integer values of this index. For instance, we obtain 

= Ki(x)--iro(x) + 0(£2), 
K2-e{x) = K2{x)--K,{x)-^Ko{x) + 0{e^) ... (60) 

X X"^ 

For the more general cases one can use the integral representation 

K,{z) = / cosh(z/t)e-"^°^^*rft (61) 
Jo 

for the index expansion. However, one can only obtain numerical results in this case. In 
expanding the integrand, we obtain 

oo ^2n 

oo ^2n oo ^2n+l 

Kr-s{z) = E(^a2„(.)-E^^^W(.) (62) 



where 



oo 

k —z cosh t , 



Jo 

POO 

aJz) = / t'' coshte-' ^°'^'dt, 
Jo 

POO 

hj,{z) = / t'' sinh 1 6-'"°'^' dt. (63) 







Due to the recurrence relations in the index u for the family Ki,{z) (cf. Sec. 15. 6^ the two formulas 
in Eq. ()62|1 suffice to calculate the Bessel function K^{z) for any u. The two functions ak{z) 
and bk{z) in Eq. (jU^ are in turn related to the functions fk{z), 

d k 
ak{z) = -—fk{z), bk{z) = -fk-i{z). (64) 
dz z 
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Using these relations, we obtain 

f2n{z). (65) 



The family of functions fk{z) is rather close to the original set of Bessel functions Ky{z) 
(cf. Fig. 13]) and can easily be studied both analytically and numerically. The limits at z ^ 
and z ^ oo are known analytically and are simple. They allow for an efficient interpolation for 
intermediate values of the argument. The functions fk{z) satisfy the differential equation 

related to the Bessel differential equation. The small z behaviour 



can be found by using yet another representation for the function fk{z) 



fk{z) = / \n\u + V^F^)du (68) 

Ji yu^ — 1 



or directly from the behaviour of the function Ky[x) at small x, 



K.^{x) = i sinh {e ln(l/x)) . (69) 
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While in four-dimensional space-time the massive propagator contains the Bessel function 
Ki^^{z), the basic function in D = 2 dimensional space-time is the Bessel function K-s{z) since 
the propagator reads 

D{x,m)\D=2 = ^Ko{mx). (70) 

As an example for the numerical evaluation of the e-expansion using Bessel functions we consider 
a toy model integral related to the one-loop case in two dimensions. We select this example 
because the integral is finite and analytically known, so that we can compare our numerical 
calculation with the exact answer. Using Eq. ()27|). we obtain 

2 J K_,{x)K_,{x)x dx = r(l + e)T{l - e). (71) 

The expansion in e is given by 

r(i - .)r(i + .) = = 1 + ^ + ^ + oie') (72) 

On the other hand, we can use the expansion 

2 4 

K.,{x) = fo{x) + '-h{x) + y,{x) + 0(6') (73) 
to rewrite the integral in the form 

poo poo 

2 / foix)'^x dx + 26"^ / fo{x)f2{x)xdx 
Jo Jo 

+77/ fo{x)U{x)xdx + - f2{xfxdx + 0{e''). (74) 
o JO 2 Jo 

Using the explicit expressions for the functions fk we have checked by numerical integration 
that the identities 

2 / foixYxdx = 1, 
Jo 

2 / fo{x)f2{x)xdx = —, 
Jo b 

1 poo 77]-4 

7: {fo{x)U{x)+3f2{xf)xdx = — (75) 
JO 3oU 

are valid numerically with very high degree of accuracy. We have implemented our algorithm 
for the e-expansion of sunrise-type diagrams as a simple code in Wolfram's Mathematica system 
for symbolic manipulations and checked its work-ability and efficiency. 

When analyzing the e-expansion one realizes that integrals over products of modified Bessel 
functions of the second kind result in integrals over products of functions fk{z). But because 



00 



2 / K_,{x)K^,{x)xdx 
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the analytical behaviour of the functions fk{z) is quite similar to the original Bessel functions, 
the numerical integration is again easy to perform. For example, in the case of the four-loop 
bubble diagram, the integral 

n4(p2 = 0) = y D{x,mfd^x (76) 
can be calculated by subtracting the expansion (cf. Appendix |E)) 

^ 1 \^,,, , ((^\' r(A) ^x^2Ar(i-A)' 



4^A+i3,2A|^ V / \^^2/ 1-A V2/ A ; 

/ /xy r(A) _ x2A r(i-A) \ 1 

^V2y l^V2y 2(1A)(2-A) 2 A(A + l)y'j ^ ^ 

from each of the propagators except for one which we keep as IR regulator. We obtain 

n4(0) = j D{x,m)(D{x,m) - A{x) + A{x)fd^x 

D{x, m) ( {D{x, m) - A{x)f + 4 {D{x, m) - A{x)f A(x) 
+6 (A(a;, m) - /\{x)f /\{xf + 4 {D{x, m) - A(x)) /\{xf + A(x)^)d^x. (78) 

The last two terms in Eq. (ff^ can be integrated analytically and can be expanded to an 
arbitrary order in e. They contain the singular contributions, i.e. all poles in e and are express- 
ible through Euler's F-functions. As expected, the pole part coincides with the expression in 
Eq. (j2H). Because the analytical expression (as given up to order in Appendix IF)) is rather 
lengthy, we present only its numerical evaluation. 



A/^nf ^(0) = - 2.5e"^ - 11.6666667e-3 - 31.701389e~2 - 67.528935e"^ 

-15871.965743 - 142923.10240^ - 701868.647626^ 

-2486982.5547^^ + 0{e^)) . (79) 



The remaining first three terms in Eq. ()78|) can be integrated numerically for D = A using the 
expansion of the Bessel functions in terms of the functions ak{z) and hk{z) |141j . No regular- 
ization is necessary since these contributions are regular at small x. The analytical expression 
for the functions to be integrated is again rather long. The zeroth order e-coefficient is found 
in Appendix O (see also the discussion of the integration procedure given in Appendix lU]) . 
However, as shown in Fig. the integrands themselves show a very smooth behaviour which 
renders the numerical integration rather simple. We obtain 

MiIiT'^{<S) = (15731. 745122 + 142349.56687e + 699112.42072^2 + 2468742.6339e^ + 0(e^)). 

(80) 
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Figure 5: Integrands for numerical integration in case of the four-loop bubble at different orders 
0(e") in € 

The sum of both parts gives 

^64(0) = m''(-2.5£-^- 11.6666667e"^- 31.701389£"2- 67.528935£"^ 

-140.220621 - 573.535535 - 2756.226905^ - 18239.9208e^ + 0(e^)) (81) 

which confirms the known result |137j . We have not included as many digits as are given in 
Ref. [TTfj . 

It is not difficult to extend the analysis to higher orders in e or to a larger number of 
significant digits in the coefficients of the e-expansion. However, since the technique is rather 
straightforward and simple we do not consider it worthwhile to extend the calculations into these 
directions. If the need arises, the potential user can tailor and optimize his or her programming 
code to obtain any desired accuracy and/or order of the e-expansion. In our evaluation we have 
used standard tools provided by Wolfram's Mathematica system for symbolic manipulations 
which allows one to reliably control the accuracy of numerical calculations. Even at this early 
stage of improvement it is obvious that our algorithm is extremely simple and efficient. 

3.8 ^-expansion of the five-loop vacuum bubble 

In this paragraph we present results for the next-order sunrise- type diagram with = 0, the 
five-loop vacuum bubble (see Fig. We have chosen to extend our calculation to the five- 
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Figure 6: A five loop vacuum bubble 

loop case since there exist no results on the five-loop bubble in the literature. The integral 
representation of the five-loop bubble with equal masses m is given by 

n5(p2 = 0) = y D{x,mfd^x. (82) 

Evaluating numerically the analytical part one obtains 

7V5fl^"''(0) = m^(3s-^ + 16.5£-^ + 51.95833£-^ + 125.6715£-2 + 259.9876£-^ 
-1360392.5934 - 16888723.177s - 111392297.46£^ 

-518606741. l£=^ + 0(e^)) (83) 

while the numerical integration of the nonsingular part gives 

^^^^^(O) = m^(l360739.9485+ 16886269.6835 

+ 111360751.916^ + 518295438.06^ + 0(£^)) . (84) 

The sum of both contributions is given by 

A/'5n5(0) = m^(3£-^ + 16.5£-^ + 51.95833£-^ + 125.6715£-^ + 259.9876£"^ 

+347.3551 - 2453.494£- 31545.55£2- 311303.l£''^ + O(£^)) . (85) 

One observes huge cancellation effects between the terms obtained by the analytical calcu- 
lation and the numerical integration. Apparently the subtraction procedure chosen here is 
non-optimal. As mentioned before, the subtraction procedure should really be optimized for 
any given problem in order to avoid a necessity to retain high numerical precision at interme- 
diate steps of the calculation. Nevertheless, our non-optimized simple subtraction procedure 
already works quite reliably with available standard computational tools. 
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Figure 7: Example of a multi-body decay of a quark in the Standard Model 



3.9 The spectral density, the discontinuity across the physical cut 
and the evaluation of the mult i- particle phase space 

Finally, we discuss the spectral density which can also be interpreted as the phase space volume 
of a multi-body decay of a single particle. An example of such a multi-body decay is shown 
in Fig. [7| At a given value of s the value of the spectral density gives the phase space volume 
of the system of particles corresponding to the lines of the sunrise diagram. The phase space 
volume for + particles with masses mi, ... , m„+i is computed according to |142j 

$„^i(p;pi, . . . = / . . . J^^^^f\Pi + ■■■+ Pn+i - P) 

where Ei = \Jpf + nil. The optical theorem relates the imaginary part of a Feynman amplitude 
to on-shell propagation of its internal lines and therefore to the product of an outer state with 
its CO nj ungate via so-called cutting rules |143j (for a recent review see also Ref. [144j). Therefore 
the discontinuity across a physical cut of the sunrise-type diagram with (n + 1) internal lines 
carrying (n+ 1) masses mi, . . . , m„+i leads to the (n+ 1) particle phase space in a natural way. 
We obtain 

*n+i(p;Pi, • • ■ ,Pn+l) = 27rp(s)|p2=^ 
which means that, up to a constant factor, the phase space is given by the spectral density. 
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4 Asymptotic analysis of sunrise-type diagrams 



In multi-loop calculations it is very useful to use expansions in different regimes of its kinemati- 
cal variables and masses, corresponding each to the smallness of one of the parameters involved 
in the problem. After defining an appropriate scale, the separation of scales can considerably 
simplify the description in the given region of energy or momentum transfer for a particular 
physical process. Such an analysis is needed when, starting from an underlying full theory, one 
wants to formulate an effective theory for a given limited energy scale. Taking for instance 
Quantum Chromodynamics (QCD) as the underlying theory, the best-known of such effective 
theories is the Heavy Quark Effective Theory (HQET) which analyzes the behaviour of a heavy 
quark near its mass shell which leads to a simplified description of physical processes subject 
to this approximation |145j . The effective theory is directly constructed through an expansion 
of the expression originally written in QCD. Its analysis, therefore, relies on the expansion of 
Feynman diagrams in a special regime. Another important example is the expansion near the 
production threshold, leading to an effective theory termed nonrelativistic QCD (NRQCD). 
NRQCD is a more complicated case of an effective theory since one introduces new variables in 
NRQCD which differ from those defined in QCD. The most involved example is Chiral Pertur- 
bation Theory (ChPT) where the variables are Goldstone modes which are not directly related 
to QCD (and where it is still difficult to construct a direct correspondence). Nevertheless, the 
common thread of all these effective theories is that they use perturbative expansions in terms 
of diagrams in given specific kinematical regimes. 

When sunrise-type diagrams enter these analyses, they usually provide the leading order 
contribution or serve as laboratory for checking more elaborate techniques. In this connec- 
tion the expansion of sunrise-type diagrams in configuration space happens to be very efficient 
for many different regimes. Having the closed form expression in Eq. ()18|1 for the correlation 
function II(p) at hand, different asymptotic regimes in the parameters (masses, external mo- 
mentum) can easily be analyzed by calculating the corresponding expansions. The first three 
subsections in this section will deal with examples in the Euclidean domain. Starting with the 
fourth subsection, we will consider some limiting cases in the Minkowskian domain which is the 
domain of physical states. The calculation of asymptotic cases is more involved in this case. 
The calculations are performed for the spectrum given by Eq. fl24j) which contains the physical 
content of the problem. 
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s-quark 




Figure 8: Leading order contribution to the two-point correlator in the sum rule analysis of 
light pentaquark states 

4.1 Large momentum expansion: close to the massless limit 

Loop integrations considerably simplify for massless internal lines. This holds true for the 
momentum space representation as well as for the configuration space representation of the 
loops. In some cases particle mass corrections to the massless approximation can become 
numerically important. For instance, in QCD the analysis of hadronic processes involving the 
strange quark require some special care since the mass of the s-quark is not small compared to 
the masses of the lightest quarks. On the other hand, the s-quark mass is still much smaller 
than the QCD scale Aqcd- The expansion in the ratio of the s-quark mass to the relevant scale 
of the considered process is therefore very effective. This expansion is especially fruitful for the 
analysis of tau lepton decays |146| I147| I148j . 

Mass corrections to the large behaviour in the Euclidean domain (i.e. expansions in mf/p^) 
are obtained by expanding the massive propagators under the integration sign in Eq. |T8|) in 
terms of the masses m,. The final integration is performed by using the identity 



Note that all these manipulations are straightforward and can be easily implemented in a system 
of symbolic computations. Some care is necessary, however, when poles of the F-function are 
encountered which reflect the presence of artificial infrared singularities. The framework for 
dealing with such problems is well-known (see e.g. Refs. |149^ IISO"] ). 

A further physics application concerns the study of pentaquark properties within the QCD 
sum rule approach. In Fig. |H1 we show the leading order contribution to the relevant two-point 
correlator of the interpolating pentaquark currents which contains one massive s-quark line. 
One therefore has to determine the mass corrections to the two-point correlator. Another 
application of the proposed techniques concerns the analysis of certain scalar mesons which can 




(86) 
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be viewed as four-quark states involving the s-quark. For instance, the scalar meson ao can be 
considered to be a KK molecule built from two i^-mesons Yet another example is the 

analysis of K*^ — mixing in the two-point function approach |152j . Note, however, that the 
approach based on the three-point function analysis appears to be more direct 

4.2 Small momentum expansion: close to the vacuum bubble 

The basic expression in Eq. (|T8|) is also well suited for finding the power series expansion in 
because the values of the polarization function and its derivatives at = can easily 
be obtained. The convenience of a p^-expansion is demonstrated by making use of the basic 
identity for differentiating the Bessel function (cf. Eq. ()A6jl ). 



Differentiation, therefore, results in an expression which has the same functional structure as 
the original function. This is convenient for numerical computations. Note that sufficiently 
high order derivatives become UV-finite. Operationally, this is obvious since the subtraction 
polynomial vanishes after sufficiently high derivatives are taken. However, this can also be seen 
explicitly from Eq. (|H7|) where high powers of suppress the singularity of the product of 
propagators at small x. The key integral is given by 

x^-^K,{Mx)dx = 2'^-2M-^r (^^) r (^^) ■ (88) 

The main physical applications are given by calculations within ChPT and HQET. The tech- 
nique is especially useful for the analysis of form factors. Another application is given by 
expansions of heavy quark correlators through the conformal mapping in order to find an ap- 
proximate spectral density from the calculation of moments |154[ll55]ll56j . Finally, applications 
can be found in the analysis of some specific rules for massive exotic states by calculating the 
moments of the spectral density. 



4.3 Dominating mass expansion: near the static limit 

Consider a hadron with one large mass which is of the same order as the external momentum 
and all other masses are small. An example is a heavy meson consisting of the c-quark and 
light quark(s). In the Euclidean domain the polarization function is not difficult to compute in 
this limit because of the simplicity of the configuration space representation and the high speed 
of convergence of the ensuing numerical procedures. One can do even better since this special 
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limit can be done analytically. When expanding the propagators in the limit of small masses 
one encounters powers of x and In(mx) (or, within the framework of dimensional regularization, 
non-integer powers of x). The remaining functions are the weight function (Bessel function of 
the first kind) and the propagator of the heavy particle with the large mass M which is given 
by the modified Bessel function of the second kind. The general structure of the terms in the 
series expansion that contribute to n(p) is given by 

^^^^^ [y) ^ •J>^(P^)^-(^^ydx (89) 
The integrations in Eq. (jH^ can be done in closed form by using the basic integral representation 

POO 

/ x''Jx(px)KJMx)dx 
Jo 

/r((A + fi + u + l)/2)r((A + fi-u + l)/2) 
2i-A'M^+A'+ir(A + 1) 
X 2Fi[{X + + U + l)/2, {X + IJ-U+ l)/2; A + 1; -p^M^) (90) 

where 2Fi{a,b; c; z) is the hypergeometric function jl35j . The corresponding integrals contain 
integer powers of logarithms. They can be computed by differentiation with respect to fi. Note 
that the maximal power of the logarithm is determined by the number of light propagators 
and does not increase with the order of the expansion: any light propagator contains only one 
power of the logarithm, as can be seen from the expansion of the Bessel function Ky{^) at small 

Physical applications in this regime are mainly sum rules for hadrons including a heavy 
flavour |157j . These are the sum rules for the Ac or A;, baryons where the lowest order contri- 
bution to the sum rules are given by the genuine two-loop sunrise-type diagram. In the case of 
the D and B mesons only the one-loop (degenerate sunrise) diagram contributes. 

4.4 Threshold expansion: close to nonrelativistic physics 

The threshold region is very important for the description of heavy quarkonia states. In the 
threshold region the necessary calculations considerably simplify if one uses an appropriate 
effective theory. This effective theory is designed to remove the degrees of freedom that can 
be treated perturbatively from the outset and considers only the dynamics of the essential 
modes relevant for the considered energy scales. For heavy quarkonia near the threshold the 
effective theory is constructed on the basis of the non-relativistic approximation for the strong 
interaction and is called nonrelativistic QCD (NRQCD). Due to the use of NRQCD there has 
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been a great advance in describing the heavy flavour production near threshold during the last 
years [T^ IT^ ITO ITTm ITT?^ ITTim ITO 

The threshold region of a sunrise-type diagram is determined by the condition + ~ 
where q is the Euclidean momentum and M = Yl,i^i is the threshold value for the spectral 
density. We introduce the Minkowskian momentum p defined by = — which allows for 
an analytic continuation to the physical cut. Operationally this analytic continuation can be 
performed by replacing q ip. To analyze the region near threshold we use the parameter 
A = M — p which can take complex values. The parameter A is more convenient in the 
Euclidean domain while the parameter E = —A = p — M is the actual energy counted from 
threshold and is used in phenomenological applications. The analytic continuation of the Fourier 
transform to the Minkowskian domain has the form 

fl{p) = 27r^+^ (^^) ^ Jx{ipx)U{x)x^^+^dx. (91) 

For the threshold expansion we have to analyze the large x behaviour of the integrand, corre- 
sponding to the region which saturates the integral in the limit p — M or, equivalently, E ^ 0. 
It is convenient to perform the analysis in a basis where the integrand has a simple large x 
behaviour. The most important part of the integrand is the Bessel function J\{ipx) which, 
however, contains both rising and falling branches at large x. This resembles the situation with 
the elementary trigonometric functions sin(z) and cos(z) to which the Bessel function J\{z) is 
rather close in a certain sense. Indeed, the function cos{z) and sin(2;) are linear combinations 
of exponentials, namely 

cos(2) = ^ (e'' + e-'') (92) 

and has also both rising and falling branches at large pure imaginary arguments: the exponen- 
tials show simple asymptotic behaviour e^^ at z = ±ioo. The analogous statement is true for 
J\{z) which can be written as a sum of two Hankel functions, 

Uz) = \{Ht{z) + H^{z)) (93) 

where H^(z) = J\{z) ± iY\{z). The Hankel functions H^(z) show a simple asymptotic be- 
haviour at infinity, 

(cf. Eq. ()A16|) ). Accordingly we split up n(p) into n(p) = n+(p) + Il~{p) with 

n±(p) = TT^+i (^^) ^ H^{tpx)U{x)x^^+^dx. (95) 
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Figure 9: Representation of Il~^{p) as a standard vacuum bubble with an additional line. The 
cross on the line denotes an arbitrary number of derivatives of that line. 

The two parts n^(p) of the polarization function Il{p) have a completely different behaviour 
near threshold which allows one to analyze them independently. This observation makes the 
subsequent analysis straightforward. 

We first consider the contribution of the n^(p) part which reduces to a regular sunrise- type 
diagram. Indeed, using the relation given by Eq. (jHH), 

K,{z) = ^e^W2^+(,^) (96) 

between Bessel functions of different kinds one can replace the Hankel function H^{ipx) with 
the Bessel function K\{px). Since the propagator of a massive particle (massive line in the 
diagram) is given by the Bessel function up to a power in x, this substitution shows that the 
weight function behaves like a propagator of an additional line with the "mass" p. The explicit 
expression is given by 

n+(p)= ^ I Ii+{x)x'^^'dx. (97) 

The function n+(a;) = Tl{x)D{x,p) is the polarization function of a new effective diagram 
which is equal to the initial polarization function multiplied by a propagator with p as mass 
parameter. We thus end up with a vacuum bubble of the sunrise-type type with one additional 
line compared to the initial diagram (see Fig. El). All derivatives of n"*'(p) = n^(M — A) with 
respect to A are represented as vacuum bubbles with one additional line carrying rising indices. 
Such diagrams can be efficiently calculated within the recurrence relation technique developed 
in |H7| EHl innj- They possess no singularities at the production threshold p = M. This can 
be seen by looking at the expansion for large x. The behaviour at large x is given by the 
asymptotic form of the functions for which one obtains 

H+(ipx) = J^e-P"(l + O(x-M), K(mx) = ./^-e-™"(l + 0(x-^)) (98) 
y znpx V 2mx 

(cf. Eqs. ^KT^ and dMH)). The large X range of the integral (above a reasonably large cutoff 
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parameter A) has the general form 

n+(M-A)~/ a;-'^e-(2*'-^)"da; (99) 

where 

a= (r2-l)(A + l/2). (100) 

The right hand side of Eq. (j^J^ is an analytic function in A in the vicinity of A = 0. It exhibits 
no cut or other singularities near the threshold and therefore does not contribute to the spectral 
density. 

In contrast to the previous case, the integrand of the second part n~(p) contains H~{ipx) 
which behaves like a rising exponential function at large x, 

(cf. Eq. ()A16|) . Therefore, the integral is represented by 

nx(M-A)~/ x'^'e-^'^dx. (102) 

The function I\~{M — A) is non-analytic near A = because for A < the integrand in 
Eq. ()l()2j) grows in the large x region and the integral diverges at the upper limit. Therefore 
the function which is determined by this integral is singular at A < (i? > 0) and requires 
an interpretation for these values of the argument A. The function is analytic in the complex 
A-plane with a cut along the negative axis. This cut corresponds to the physical positive energy 
cut. The discontinuity across the cut gives rise to the non- vanishing spectral density of the 
diagram (cf. Appendix . 

The integral for T\~{p) in Eq. ()95|) cannot be done analytically. In order to obtain an 
expansion for the spectral density near the threshold in analytical form we make use of the 
asymptotic series expansion for the function 11 (x) which crucially simplifies the integrands but 
still preserves the singular structure of the integral in terms of the variable A. The asymptotic 
series expansion of the main part of each propagator, i.e. of the modified Bessel function of the 
second kind, to the order is given by Eq. ()A15jl . The asymptotic expansion of the function 
n(x) therefore consists of an exponential factor e~*^^ and an inverse power series in x up to 
the order A^, where A^ is closely related to A^. It is this asymptotic expansion that determines 
the singularity structure of the integral. We write the whole integral in the form of a sum of 
two terms, 

n-(p) = TT^+iy (^^^ ^iJ-(ipx)(n(a;) -n^^(x))x2^+^+2^da; (103) 
+7r^+i / (^) ' ^A"(w)n^'(x)x'^+^+'^rfx = t[''\p) + ti'^^p). 
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The integrand of the first term 11'^* (p) behaves as at large x while the integrand of the 

second term accumulates all lower powers of the large x expansion. Note that only the large 
X behaviour is essential for the near threshold expansion of the spectral density. This fact has 
been taken into account in Eqs. (jHSl) and ()l()2j) where we introduced a cutoff A. However, from 
the practical point of view the calculation of the regularized integrals with an explicit cutoff is 
inconvenient. The final result of the calculation - the spectral density of the diagram - is inde- 
pendent of the cutoff, but the integration is technically complicated if the cutoff is introduced. 
However, in extending the integration over the whole region of the variable x without using 
the cutoff one immediately encounters divergences at small x because the asymptotic expan- 
sion is invalid in the region near the origin, so one is not allowed to continue it to this region. 
The standard way to cope with such a situation is to introduce dimensional regularization. It 
allows one to deal with divergent expressions at intermediate stages of the calculation and is 
technically simple because it does not introduce any cutoff and therefore does not modify the 
integration region very much. Note that dimensional regularization does not necessarily regu- 
larize all divergences in this case (in contrast to the standard case of ultraviolet divergences) 
but it nevertheless suffices for our purposes. We shall use a parameter e to regularize the 
divergences at small x. 

The first part H°'*(p) in Eq. containing a difference of the polarization function and 

its asymptotic expansion since the integrand gives no contributions to the spectral density up 
to a given order of the expansion in A. This is so because the subtracted asymptotic series to 
order cancels the inverse power behaviour of the integrand to this degree N . The integrand 
decreases sufficiently fast for large values of x and the integral converges even at A = 0. 
Therefore this term is inessential for the evaluation of the expansion of the spectral density up 
to given order. 

The expansion of the spectral density at small E is determined only by the integral H"'*(p) 
in Eq. ()103j) . This integral is still rather complicated to compute but we can go a step further 
in its analytical evaluation. Indeed, since the singular behaviour of H"*(p) is determined by 
the behaviour at large x, we can replace the first factor, i.e. the Hankel function, in the large 
X region by its asymptotic expansion up to some order N . We use Eq. (IA16|1 to obtain a 
representation 

II^-(p) = 7r^+i|(^!^^ iJ--(^px)H^^(x)a;2^+^+2"rfa;. (104) 

The index "rfas" stands for "double asymptotic" and indicates that the integrand in Eq. ()l()4j) 
consists of a product of two asymptotic expansions: one for the polarization function H(x) 
and another for the Hankel function H\{x) as weight (or kernel). Both asymptotic expansions 
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are straightforward and can be obtained from standard handbooks on Bessel functions (cf. 
Appendix El- We therefore arrive at our final result: the integration necessary for evaluating 
the near threshold expansion of the sunrise-type diagrams reduces to integrals of Euler's Gamma 
function type, i.e. integrals containing exponentials and powers. Indeed, the result of the 
expansion in Eq. ()104|) is an exponential function e~^^ times a power series in l/x, namely 

where a has already been defined in Eq. and the coefficients Aj are simple functions of the 

momentum p and the masses mj. The expression in Eq. ()105|) can be integrated analytically 
using 

/•OO 

/ x-"+'"e-^"rfx = r(l - a + 2£)A"-^-2^ (106) 

The result is 

Af-l 

j^das^^ _ A) = ^ A^Y{\ -a-j + 2£)A"+^-i-2e_ ^^Q7) 

3=0 

This expression is our final representation for the part of the polarization function of a sunrise- 
type diagram necessary for the calculation of the spectral density near the production threshold 
(see Appendix The spectral density is a function oi E = —A and will be denoted by 
p{E) = p{{M + E)"^) in the following. Starting from the main result in Eq. (jl(J7|) . we discuss 
the general structure in detail. In the case where a takes integer values, these coefficients result 
in divergences for small values of e. The powers of A in Eq. (|1U7|) have to be expanded to 
first order in e and give 

— A^^ = — + lnA + 0(£). (108) 
2£ 2e 

Because of 

Disc ln(A) = \n{-E - iO) - \n{-E + iO) = -2m9{E), (109) 

n'^"*(M — A) in Eq. p07p contributes to the spectral density. For half-integer values of a the 
power of A itself has a cut even for 6 = 0. The discontinuity is then given by 

Disc TA = -2zv^ 9{E). (110) 

Our method to construct a threshold expansion thus simply reduces to the analytical calculation 
of the integral in Eq. ()104|) which can be done for arbitrary dimension and an arbitrary number 
of lines with different masses. In the next paragraphs we use our technique to work out some 
specific examples which demonstrate both the simplicity and efficiency of our method. 
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Figure 10: Various results for the spectral density for n = 3 equal masses in Z) = 4 space- 
time dimensions as function of the threshold parameter E/M. Shown are the exact solution 
obtained by using Eq. ()112j) (solid curve) and threshold expansions for different orders taken 
from Eq. ()114|) (dashed to dotted curves). 



4.5 Equal mass genuine sunrise diagram at threshold 

The polarization function corresponding to the genuine sunrise diagram with three propagators 
with equal masses m in D = 4 space-time dimensions is given by 



The exact spectral density is given by the integral representation in Eq. (j24jl which for this 
particular case reads 

p{s) = / h{xy/s)U{x)x^dx. (112) 

l\/ S J c—ioo 



In order to obtain a threshold expansion of the spectral density in Eq. ()112|) we use Eq. pOTj) 
to calculate the expansion of the appropriate part of the polarization function. To illustrate 
the procedure we derive the explicit functional behaviour of the integrand in Eq. p04|) which 
is given by an asymptotic expansion at large x, 

X |l+ ^ 3^ 9 27 15 -3x1 /^^gx 

1 8mx 8px 128m^a;^ GAmpx"^ 128p^a;^ J 
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From Eq. ()113|) we can easily read off the coefficients Aj that enter the expansion in Eq. ()105|) . 
The spectral density is obtained by performing a term-by-term integration of the series in 
Eq. ()113|) and by evaluating the discontinuity across the cut along the positive energy axis 
E > 0. The result reads 

p(E) = ^{1 ?7H rf - -rf ^ ri* ri^ (114) 

^ ' 3847r3v^\ 2 ' 16 ' 8 ' 128 ' 256 ' ^ ' 

129 . 3 7 4047 o 18603 „ 248829 ^ 

- ~ ^^^^^^ + T^^V - r-o.ooo ^ + 0{v ) \ 



1024 ' 256 ' 32768 ' 65536 ' 524288 

where the notation rj = E/M, M = 3m is used. The simplicity of the derivation is striking. 
By no cost it can be generalized to any number of lines, arbitrary masses, and any space-time 
dimension. The genuine equal mass sunrise has been chosen for definiteness only. It also allows 
us to compare our results with results available in the literature. Eq. ()114|) reproduces the 
expansion coefficients dj obtained in Ref. |168j (the fourth column in Table 1 of Ref. |168j ) by 
a direct integration in momentum space within the technique of region separation |169j . In 
Fig. El the exact solution is shown together with expansions to various orders. 

4.6 Three loops and a route to any number of loops at threshold 

The sunrise-type diagrams with four or more propagators cannot be easily done by using the 
momentum space technique because it requires the multi-loop integration of entangled mo- 
menta. As emphasized before the configuration space technique allows one to immediately 
generalize our previous results to any number of lines (or loops) with no additional effort. Con- 
sider first the three-loop case of sunrise-type diagrams. The polarization function of the equal 
mass sunrise-type diagram with four propagators in D = 4 space-time is given by 

nw^^^. (115, 

The exact spectral density of this diagram can be obtained from Eq. (jll2j) while the near 
threshold expansion can be found using Eq. p07j) . The expansion of the spectral density near 
threshold reads 

E'^'^M^/^ 1 81 2 2811 3 17581 4 

piE) = r7H f]^ r7'^H (116) 

268807r5v^\ 4 ' 352 ' 18304 ' 292864 ' ^ ^ 

1085791 , 597243189 « 4581732455 7 496039631453 « 

H ri^ ri^ H ri^ ri^ + Oirf) \ 

19914752' 3027042304' 12108169216' 810146594816' ^''] 

where rj = E/M and M = Am is the threshold value. One sees the difference with the previous 
three-line case. In Eq. ()116j) the cut corresponds to a square root branch cut while in the 
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three-line case one has a logarithmic cut. One can easily figure out the reason for this by 
looking at the asymptotic structure of the integrand. For an even number of lines (i.e. odd 
number of loops) one has a square root branch cut, while for an odd number of lines (even 
number of loops) one has a logarithmic branch cut. This statement generalizes from D = 4 to 
any even space-time dimension. In the general case the structure of the cut depends on the 
dimensionality of the space-time. The general formula for the ra-loop case reads 

For D = 4 space-time dimension (i.e. A = 1) we can verify the result of Ref. jHZlEHlEni, 

p(E) ~E(3"-2)/2(i + o(E)). (118) 

Returning to Eq. ()116|1 one has numerically 

p{E) = 8.5962- 10"^E^/V/^{l.000 -0.250r/ + 0.230r72 (119) 
-0.15V + O.OeOr;^ + 0.055?]^ - 0.197//*^ + 0.3787^^ - 0.612r/^ + 0{ri^)} 

where we have written down the coefficients up to three decimal places. It is difficult to say 
anything definite about the convergence of this series. By construction it is an asymptotic 
series. However, we stress that the practical (or explicit) convergence can always be checked 
by comparing series expansions like the one shown in Eq. ()119|) with the exact spectral density 
given in Eq. ()112|1 by numerical integration. 

We conclude that the spectral density of the sunrise-type diagram can be efficiently calcu- 
lated within the configuration space technique. It does not matter whether one is aiming for 
the exact result or an expansion the configuration space technique which can readily deliver 
the desired result. The exact formula in Eq. ()112|1 as well as the threshold expansion obtained 
from it can be used to calculate the spectral density for an arbitrarily large number of internal 
lines. We stress that the case of different masses does not lead to any complications within 
the configuration space technique: the exact formula in Eq. ()24j) and/or the near threshold ex- 
pansion work equally well for any configuration of masses. We do not present plots for general 
cases of different masses because they are not very illuminating as one can only see the common 
threshold. However, there is an interesting kinematic regime for the different mass case which 
is important for certain applications which, to the best of our knowledge, has not been treated 
in the literature before. An analytical solution for the expansion of the spectral density in this 
regime is given in the next subsection. 
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4.7 New features of the threshold expansion: resummation of small 
mass effects for strongly asymmetric mass configurations 

The main question in constructing expansions is their practical usefulness and the region of 
validity. The threshold expansion for equal or almost equal masses breaks down for instance 
for E M = ^nii. However, if the masses are not equal, the region of the breakdown of the 
expansion is determined by the mass with the smallest numerical value. The simplest example 
where one can see this phenomenon is the analytical expression for the spectral density of the 
simple loop (degenerate sunrise- type diagram) with two different masses mi and m2- In D = 4 
space-time dimensions (see e.g. Ref. jHHEHlEni) one has 



lE{E + 2m,){E + 2m2){E + 2M) 

= {Am+E)f 

where M = mi + 1112. The threshold expansion is obtained by expanding the right hand side of 
Eq. ()120|) in E for small values of E. If 1712 is much smaller than mi, the expansion breaks down 
at ~ 2m2. The break-down of the series expansion can also be observed in more general 
cases. If one of the masses (which we call mo) is much smaller than the other masses, the 
threshold expansion is only valid in a very limited region E < 2mo. 

To generalize the expansion and extend it to the region of -E ~ M one has to treat the 
smallest mass exactly. In this case one can use a method which we call the resummation of 
small mass effects ^O]- This method will be explained in the following. We start with the 
representation 

flP-(p) = TT^+^I 'i7,--(zpx)n- (a;)^2^+^+^^da; (121) 

which is the part of the polarization function contributing to the spectral density. The integrand 
in Eq. ^\2\\) has the form 

n-(x) = n-i(a;)D(mo,a;) (122) 

where the asymptotic expansions are substituted for all the propagators except for the one with 
the small mass mo. This is indicated by the index ";?as" in Eq. (jl21|) which stands for "partial 
asymptotic". The main technical observation leading to the generalization of the expansion 
method is that W°''^{p) is still analytically computable in a closed form. Indeed, the genuine 
integral to compute has the form 
00 

v/^(2/3)-r(/x + z/)r(/i-z.) ^ ffx + u fx + u + l 1 p'\ 
- (25)^+- r(/i + l/2) ^^^1,^' 2 '^+2'^"^j 
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where a = A — mo and (3 = mo. The integral W^^lp) in Eq. p21|) is thus expressible in 
terms of hypergeometric functions |17()j . For the construction of the spectral density one has 
to find the discontinuity of the right hand side of Eq. ()123|) . There are several ways to do 
this. We proceed by applying the discontinuity operation to the integrand of the integral 
representation of the hypergeometric function. The resulting integrals are calculated again in 
terms of hypergeometric functions. Indeed, 



1 

Disc / x''-^e''''KJpx)dx 

2m Jo ^ ' 



2^{a^ - pyl^-^ r(3/2) ^ ( l-^^-u 2-^,-u ^ 

ay''- 13- r(3/2-/x) 1^ 2 ' 2 ' 2 ^' j U^^^ 

where a = E + mo- The final expression in Eq. ()124j) completely solves the problem of the 
generalization of the near threshold expansion technique. For integer values of /i there are 
no singular Gamma functions (with negative integer argument). Therefore we abolish the 
regularization and set e = when using this expression. We have thus found a direct transition 
from the polarization function (expressed through the integral) to the spectral density in terms 
of one hypergeometric function for each genuine integral. There is no need to use the recurrence 
relations available for hypergeometric functions. 

As an example we present the (over) simplified case of the two-line diagram with masses m 
and mo -C m in four space-time dimensions. We cite this example because the expansion of 
the spectral density and its generalized expansion can be readily compared analytically with 
the exact result in Eq. ()120|) . Both expansions, the pure and the resummed expansion, will be 
compared with the exact anaytical result. 

The pure expansion of the spectral density near threshold (the second order asymptotic 
expansion should suffice to show the general features in a short and concise form) is given by 



-das 



^ V2^^( fl 1_7_\E 
^ ' SttHP/' \ Vm mo Mj 4 



where M = m + mo- As mentioned above, this series breaks down for E > 2mo (see Eq. p2()|l ). 
The analytical expression for the spectral density of the polarization function in Eq. ()12H) 
within the generalized asymptotic expansion based on Eq. ()124|1 is given by 



JmEiE + 2mo) j / 1.3 mg \ 

87r2(E + M)3/2 \2 '2' 2' {E + mo)y 



E{E + 2mo) /I 5 mg \ 

+ 8m(E + M) ^' 2 ' ^ " {E + mof ) ^^^^^ 
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Figure 11: Various solutions for the one- loop spectral density for two different masses m and 
mo -C m in D = 4 space-time dimensions. Shown are the exact solution which is obtained by 
using Eq. ()112|1 (solid curve), the pure threshold expansions using Eq. p25|l (dotted curves), 
and the solutions for the resummation of small mass effects such as in Eq. ()126|) (dashed curves). 
For both case the expansions are shown from the first up to the fourth order in the asymptotic 
expansion. For the pure threshold expansion the order is indicated explicitly. 



E\E + 2moY / 16m{E + M) \ 3 7 

128m^{E + My[ ^ 5(E + mo)2 ) ' '1^'2'2' (E + mo^ 

We have set the regularization parameter e = because the spectral density is finite. With 
e = the resulting expressions for the hypergeometric functions in Eq. ()124|) simplify. The first 
term in the curly brackets of Eq. (jl26p is obviously equal to 1 in this limit because the first 
parameter of the hypergeometric function vanishes for e = 0. However, we keep Eq. p26|) in 
its given form to show the structure of the contributions. The generalized threshold expansion 
has the form 

pP'^'iE) = go{E, mo) + Eg,{E, mo) + E^g^iE, mo) + . . . (127) 

where the functions gj{E, mo) represent effects of the resummation of small mass effects and are 
not polynomials in the threshold parameter E. In the simple one-loop case the hypergeometric 
functions reduce to elementary functions. For instance. 
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3(E + mo) ml , E + mo + JE{E + 2mo) 

' E + mo m ' 



2E{E + 2mo) y 2Je{E + 2mo) \E + mo - J E{E + 2mo) 



Higher order contributions are given by hypergeometric functions with larger numerical values 
of the parameters. They can be simplified by using Gaussian recurrence relations for hyperge- 
ometric functions.^ 

In Fig. ^2 we compare the exact result of Eq. ()112j) with the pure expansion in Eq. p25p 
and the resummed expansion in Eq. (jl26|) . While the pure expansion breaks down already 
at the order E ~ mo, the convergence of the expansion in Eq. p26|) breaks down only at 
E M = m + mQ. The resummation leads to an essential improvement of the convergence in 
comparison with the pure threshold expansion. Further examples and their discussion can be 
found in Ref . [H! . 

4.8 Large s expansion: Expansion of the spectrum at large energies 

In the Minkowskian domain where s = the situation is closely related to the large p expansion 
if all masses are small compared to the external momentum. We start with a simple example 
which will be calculated in momentum space as well as in configuration space. In the next 
paragraph we give more involved examples. 

As a first example for an expansion in m^/ s we look at a case where the spectral density is 
given by a finite series only. For the one-loop sunrise-type diagram with one massive and one 
massless line we can use Feynman parametrization and obtain in (Euclidean) momentum space 

fl(^p) = JM_. f\l - (xp'^ + m'^y dx. (129) 
(47r)^ ^ Jo ^ ^ 

In order to determine the spectral desity, the discontinuity can be calculated already on the 
level of the integrand. A discontinuity appears if xp'^ + m^ < which is satisfied for p"^ < — m^. 
We refer to Appendix O for more details. If we parametrize p"^ = se^'^, we can reach the cut 
from both sides cp = —n and = +tt and can calculate the difference, i.e. the discontinuity 

27ri 

Disc (m^ — xs)~^ = r(xs — m^)^^6'(xs — m^). (130) 

i {£)l yi- £) 

For the spectral density we obtain 

,(,) = J-Kscnip)l___^ = (4.)4(i-ri CJ' - - "T' (131) 



^Note that Eq. 1)1281) does not lead to the exact function in Eq. I|12U|) because terms of order are 
missing which originate from the difference part Il'^^{p) of the correlator. The difference part simply corrects 
the behaviour of the coefficient functions by the small mass contributions. 
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Because the general factor is no longer singular, we can calculate the integral for e = where 
the integrand is simply unity. Then we obtain 

This result can also be obtained by using configuration space techniques. At first sight one may 
think that this example is too simple to be worth to be calculated using configuration space 
techniques. Later on we shall see that for the multi-loop case the configuration space method 
is the only one which leads to an analytical result. 

In configuration space the spectral density of the one-loop diagram with one massive and 
one massless line in four-dimensional spacetime is given by Eq. pi2j) . The integral contains 
the modified Bessel function of the first kind, Ii{z). Because Ki{x) is a pure exponentially 
decreasing function in its asymptotic expansion while Ki{x) ± i7r/(x) increases exponentially 
along the real axis if continued in the upper resp. lower complex plane (note the ambiguity due 
to Stokes' phenomenon mentioned in Appendix [CI), we can divide according to 



in j h{z)f{z)dz = j {inh{z) - Ki{z))f{z)dz + J Ki{z)dz, (133) 
in r Ii{z)f{z)dz = r (inhiz) + Ki{z))f{z)dz - r Ki{z)dz (134) 



and can lower the paths to the negative resp. positive real axis by using Cauchy's theorem 
realizing that the quarter circle integrals in the corresponding quadrant vanish (see Ref. jSH|)- 
The origin and the negative real axis is circumvented by two half-circle paths C_ (lower half 
plane) and C+ (upper half plane) in positive direction. The result found in Ref. |SH] reads 

pc+ioo POO , , 

in \ h{x^s)Y\{x)x^dx = \ K^{r^s)[m{r) -J\{e^r) -J\{e-'^r)\r^dr (135) 

+ f {inli{zy^) + Kl{z^/s))U{z)z^dz+ f {inIl{z^/^) - Kl{z^/s))U{z)z^dz. 
JC- Jc+ 

For our simple example we have to take the correlator function n(a;) = D{x, m)D{x, 0). While 

the massless propagator is invariant under the multiplication of e^*'^ = ±1 of its argument, i.e. 

D{—x, 0) = D{x, 0), for the massive propagator we obtain 

Die^'^r.m) = / (Ki(mr) ± inli{mr)) (136) 
{2n)'^r'' 

where Kx{e^''^x) = e'^'^'^Kxix) =F inhix) is used. We obtain 2n(r) - n(e^''r) - U{e-''^r) = 

which results in the vanishing of the first part in Eq. p35|) . We are therefore left with the two 

semicircle integrals which we can combine in a different manner 

—2n7 r 2 r 

p{s) = ^^ h{z^)T{{z)z^dz- ^ Ki{z^s)T{{z)z^dz. (137) 

VS JC-+C+ y/s JC--C+ 
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In expanding the first integrand for small values of r = \z\ we obtain 



h{z^s)Ii{z)'. 



+ 0{z) 



(138) 



the integral results in 



In r + 2^9 



(139) 



For the second integrand the expansion results in many more terms 

1 1 



Ki{z^)U{z)z^ 



+ 



{27T)^z^y/s 4(27r)4;zyi 



m ln(m ) + sins) 



m + s 



2(2vr)%yi 0"©+^^"^ 



(140) 



However one has 
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■-r 
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■ -1 ■ 


c+ z^ 


UzA 


C--C+ 









r dz 
Jc--c+ z 



\nz 



C--C^ 



In e + i(f 



ip=—'K 



(/3=0 



The only non-vanishing contribution is given by 



dz Ir 2 

— in 2; = - in z 

C--C+ z 2L 



C--C+ 



IT 



7r-7r = 0. (141) 



(142) 



Therefore, we obtain 



171 ~h 5 



The spectral density, finally, reads again 

— 27ri i^/s 2 rn? + s 1 



167r3 327r2^ IGtt^ 



m 
s 



(143) 



(144) 



as in Eq. 

4.9 Large s expansion for massive one- and multiloop diagrams 

A second example that we treat is the one-loop diagram with two different masses mi and 777,2. 
We start again with the calculation in momentum space to obtain 

r(2-D/2) n 



n(p) 



(47r)^/2 



{x{l — x)p^ + (1 — x)ml + xmfj 



2\ D/2-2 



dx. 



(145) 
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Again the discontinuity can be calculated already on the level of integrands. Using Eq. ()C15|) 
we obtain 



Disc x(l — x)s + (1 — x)m^ + xml^ = (146) 
27ri 

rr 



^ (^x{l — x)s — (1 — x)m1 — xml^ 6 (^x{l — x)s — (1 — x)m\ — xm^ 



In order to write down the spectral density we have to determine the limits for x given by 
the positiveness of the theta function. The zeros of the argument of the theta function can be 
calculated to be 

1 / m? ml\ 1 r— — r— 

xi,2 = - 1 + ^ - ^ ± - VA(1, ml/s, mils) (147) 



where 

A(x, y, z) = x"^ + + - 2xy - 2xz - 2yz (148) 

is Kallen's lambda function. In terms of these two zeros xi < X2 we can write the condition for 
the argument of the theta function as 

{x — Xi){x2 — x)s > =^ Xi < X < X2. (149) 

The two zeros take real values if 

^A(l, mj/s, mj/s) = {s - (mi - ms)^) (s - (mi + ms)^) > 0. (150) 



This is the case for s > (mi + 1712)^ ot s < (mi — 1712)"^. We do not consider the part of the 
spectral density up to the pseudothreshold (mi — 1712)^ because in this case both xi and X2 are 
negative. Starting from the threshold (mi + m2)^, both xi and X2 take values between and 
1. We obtain 

1 ~ 1 p2 _ _ _ 

P^'^ = 2^°^^^n(-.) = (4,).-.r(i_,) X^ s-^ix-x,rix2-xrdx 



A(l,m?/s,mi/s) (151) 



(47r)2-^r(2 - 2e) 

where we have used the substitution x' = {x — Xi)/{x2 — Xi) and Euler's Beta function 

i?(a,;5) = £l^llM= [\-\l-xY-'dx. (152) 



r(a + p) Jo 

For D = 4, i.e. £ = 0, we obtain 



Pis) 



{An)^\ 



'2 2\ / 2 2\2 

mf mo \ / mf m. 



1-2^ + ^ +^ ^ . (153) 
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For the case m2 = we recover the spectral density of the previous example. If we consider an 
expansion in terms of large values of s compared to the threshold value s = (mi + 1712)^ , i.e. if 
we expand in small values of 2; for ^ = (mi + m2)^/s, we obtain 

The same result shall now be calculated in a different way. Consider the correlator function 
n(p) = 27r^+^ (^y) ^ Jx{px)D{x, mi)D{x, m2)x^^^^dx. (155) 

We can expand the propagators D{x, irti) in nxij \fs^ leading to pure powers in x. Starting from 
Eq. ()8(i|l . the integration of the Bessel function of the first kind with (non-integer) powers of x 
can be cast into a more appropriate form 

"-(f)"^^^(-)^-'--(l)"""^"^^^^^- (-) 

We can then immediately calculate the discontinuity and obtain 

2^1°- (y)" ^^^^'^'^ - 

Using the expansion of the Bessel functions K\(^z) at small values of 2;, the expansion of the 
propagator for A close to A = 1 for small mass m (compared to y^) reads 

^, . (mxy'K\(mx^ 
D[x,m] = 



(27r)^+ix2^ 

2^-ir(A) / 1 fmx\^ r(l-A)/mx\2^ ^ / , ,4, ,2+2x\\ /.ton 

Eq. p55|) contains two of these propagators. We now expand up to the power m^ and keep A 
close to 1. We obtain a zeroth order contribution, two contributions of the order m^, and two 
contributions of the order m?^^~^^\ We consider these contributions in turn. 



zeroth order (// = — 4A): 



27r^+i(s/4)^-i /2^-ir(A)V 271^ / r(l 



2r(2A)r(A) \{27t)^+^ J 2r(2)r(l) V(27r)2y 167r2 



(159) 



second order (/i = 2 — 4A, z = 1, 2): 

27r^+^(./4)^-^ / 2^-^r(A) y (mf/4) . . 

2r(2A-i)r(A-i) [1^) (r(A) = (A-i)r(A-i)) (leo) 

-27r^+i(s/4)^-2 /2^-ir(A)y m2 ^ -2tt\s/A)-^ / T{1) ^ mf _ -m^ 

2r(2A- l)r(A) V (27r)^+i J 1 ' 2r(l)r(l) l^(27r)^+ij ~T ~ Wn^s 
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order m^^^^^^ (also called "second order primed", fx = — 2A, i = 1,2): 

_ js/^)-' / 2^-^r(A) y 2r(l^A) (mT 

r(-A - ^/2)T{-^/2) [ (27r)A+i ) r(l + A) ^ ^ 



(161) 



Since the third contribution vanishes for /j, — > — 2A, only the first two contributions have to be 
taken into account. One ends up with 

in agreement with the result obtained earlier. Moreover, we can predict the result for a n-loop 
sunrise-type diagram with different masses mi, . . .m„+i in four space-time dimensions where 
II — —2{n + 1)A (zeroth order contribution), n — 2 — 2{n + 1)A (second order contribution), 
and // = —2n\ (second order primed contribution which does not vanish for n > 1), 

/ 2^-^r(A) \ {s/Ar~' ( 
P^'^ ^ IW^J 2r((n+l)A)r(nA) 1,^ + ^'^^^ 

((n+l)A-l)(nA-l) 1^ rr^ _ r((n + l)A)r(l - A) K V , ^ ( ]_ 
1-A ^ r{{n-l)X)r{l + X)^A^ J 

The second and third term have to be expanded in £ = 1 — A in order that the singularity in e 
cancels. After a few simplifications using expecially 

ijj{n + 1) ^ 7jj{n) + - (164) 



for the polygamma function, we finally obtain for A = 1 

^<^) - (mn-iy. (' + " I' ((" - » (t) + ^('^C') + ^-)) - ") V + <^ (>, 

(165) 

This result is also valid forn = 1 and again confirms our previous results. Note that {ip{n) +^e) 
is a rational number because Euler's constant cancels out. Indeed, we have 

k=l 

Note, however, that forn > 1 logarithmic contributions ln(mf/s) appear in the spectral density. 
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5 Non-standard propagators and other exotic settings 



In this section we deal with modifications of the standard propagators in sunrise-type diagrams 
in any space-time dimension. First we consider odd-dimensional space-time which in most of 
the cases allows one to calculate the sunrise-type integrals in closed form. Returning to space- 
time dimensions close to four, we then deal with larger powers of propagators leading to larger 
values for the Bessel function indices. In this case one can apply recurrence relations to reduce 
integrals containing higher powers of the propagators to a set of master integrals [211 • The 
master integrals themselves show Laplace-type asymptotics (for a definition see e.g. Ref. |171j ). 
Finally we deal with nontrivial numerators in two examples involving vacuum bubbles. 



5.1 Odd-dimensional case 

Compared to four-dimensional space-time, integrals containing products of Bessel functions 
in odd dimensional space-time are analytically solvable because the Bessel functions simplify 
significantly. It is interesting to note that the evaluation of Eq. (fTHj) can be done in a closed 
form for any number of internal lines in odd- dimensional space-time. One might suspect that 
this special case is too simple to have any practical applications. However, at the end of 
this subsection we will refer to many applications from various fields of physics involving odd- 
dimensional space-time. As the simplest example we take three-dimensional space-time D — > 
Do = 3. For Do = 3 (and Ao = (Do - 2)/2 = 1/2) the propagator in Eq. © reads 

D(x, m) Dsix, m) = (,^^3/2^ = (167) 

With A = Ao = 1/2 the weight function becomes 

, . 2 smjpx) 



after angular integration. The explicit result for the {n + l)-line sunrise diagram is then given 
by the integral 



n(p) = Att f '-^I^^-^i^x^dx 
^ ' Jo px (47rx)"-i^ ' 



r(l-n + 2e) 



(M - ip)''-'-'' - (M + ip)""-'-'' fi'' (169) 



2iip(47r)" 

where e is used for regularization and M = X)'^*- Note that we have used an unorthodox 
regularization method by multiplying a factor (/ix)^^. We will return to this point later on. 
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We consider some particular cases of Eq. ()169|) for different values of n. For n = we simply 
recover the (Euclidean) propagator function n(p) = (M^ + p^)~^ with the discontinuity 

= Hi^IlM = 1_ { n(p)| fl(p) ^ = - m^) (170) 

' — ' — ' p^=s exp(— t7r) / 



p2=s exp{i7r) 



where s is the squared energy, s = —p^. As remarked on earlier it is appropriate to call this 
expression the spectral density associated with the diagram. For n = 1 the answer for the 
polarization function n(p) is still finite (no regularization is required) and is given by 

~, . I , (M + ip\ , , 

Hip) = In ^ . 171 

Snip \M-ipJ ^ ' 

The spectral density, i.e. the discontinuity of Eq. (jl71|) across the cut in the complex p^-plane 
is given by 

p{s) = -^eis-M'), s = -p\ s>0 (172) 

which is nothing but the three-dimensional two-particle phase space. This can be immediately 
checked by direct computation. The cases with n > 1 have more structure and therefore are 
more interesting. For the genuine sunrise diagram with n = 2, Eq. ()169p leads to 

1 fM + ip\ J (M'^ + p^ 

327r2 U ~ ^ \M -ip J ~ ^ I ^ 



n(p) = ^(--^>.>(^l-i"(^ . (173) 



While in Eq. fll7H|l the arbitrary scale iJ? appears due to regularization, the spectral density 

is independent of this scale. This is because the spectral density is again finite and, therefore, 
independent of the regularization used. The general formula for the spectral density for any 
n > in D = 3 can be extracted from Eq. p69|) . It reads 

p(,) = /^7^C^ ^(^ - ^')- (175) 
'^^ ' 2(47r)"(r2- l)!v^ ^ ' ^ ' 

We now want to comment on the relation between the momentum subtraction and our 
unorthodox dimensional regularization. Taking Eq. ()169|) for n = 2 with momentum subtraction 
at the origin, one obtains 

n(p) = r ( !^ _ -^(^:,)^^dx (176) 



px J {Airyx 

which is UV-finite even for e = because there is no singularity at the origin. For practical 
computations it is convenient to keep the factor (fix)'^'^ in the integrand since this factor gives 
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a meaning to each of the two terms in the round brackets in Eq. ()176p . Using this factor we 
obtain 

^ ' Jo \ px J {Ait)^x^ ' 



2e 



1 f£,„f^u.„f^^;±i^n. (177) 



327r2 \ ip \M-ip) \ 

The poles cancel in this expression and the arbitrary scale /i changes to M. This corresponds 
to a transition from MS-type renormalization schemes to the momentum subtraction scheme 
(in this particular case with a subtraction at the origin). Since the spectral density p{s) is finite 
it can be computed using any regularization scheme as can be seen from Eqs. ()173|1 and ()177|1 . 

We mention that in the three-dimensional case the spectral density p{s) can also be found 
for general values of n by traditional methods since the three-dimensional case is sufficiently 
simple. By using the convolution for the evaluation of spectral densities one stays in the same 
class of functions, i.e. polynomials in the variable y/s divided by y/s. The general form of the 
convolution equation in Z)-dimensional space-time reads 

^n{s) = J ^kisi)%is2)^iis,Si,S2)dsids2, k + p + I = u. (178) 

For the particular case of three-dimensional space-time the kernel $i(p^, mf, is given by 

(27r)2$i(p2, ml ml) = J 6{k^ - mi)8{{p - kf - ml)d^k (179) 

or explicitly by 

$i(s,si,s2) = -^e{s-i,/ri + ^f). (180) 

Eq. p8()|l can be seen to be the two-particle phase space in three dimensions (cf. Eq. ()172|l ). 
This is a rather simple example. However, our technique retains its efficiency for large n. 

We list some potential applications of the general results obtained in this paragraph for odd- 
dimensional space-time. In three space-time dimensions our results can be used to compute 
phase space integrals for particles in jets where the momentum along the direction of the jet 
is fixed |lU4j . Another application can be found in phase transitions, for instance the three- 
dimensional QCD which emerges as the high temperature limit of the ordinary theory of strong 
interactions for the quark-gluon plasma (see e.g. [5^ 1^ ri72[ I92j). Three-dimensional models 
are also used to study the question of dynamical mass generation and the infrared structure of 
the models of Quantum Field Theory in general |H7 |ll73|ll74j and some problems of QCD |175j . 
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Note that particular models of different space-time dimensions are very useful because their 
properties may be simpler and may thus allow one to study general features of the underlying 
field theory. For example, in six- dimensional space-time the simplest model of quantum field 
theory 0^ is asymptotically free and can be used for simulations of some features of QCD [ 17(ij . 

As we have already mentioned before, closed form analytical results can be obtained for 
odd dimensional space-time D = 3, 5, 7, . . .. For example, the propagator in five-dimensional 
space-time (A = 3/2) reads 

D{x,m) D,ix,m) = ^"^^^^^^^^^^^^^ = £^(1 + ^^) (181) 

which assures that the basic integration can be performed in terms of elementary functions 
(powers and logarithms) again. 

Applications for odd space-time dimensions other than three can be found in some models 
in unified field theories, or also in a general analysis of the divergence structure in QFT. Five- 
dimensional models of QFT are rather popular for general purposes |177^ I178t I179| IT^ . They 
have useful applications for Yang-Mills theories in five-dimensional space-time where the UV 
structure of the models can be analyzed |181j . 

Finally, non-standard space-time dimensions are useful in order to obtain estimates for the 
standard case Dq = 4. The feature that K^{x) > K^{x) holds for modified Bessel functions 
of the second kind when v > ji and for positive x-arguments leads to K^i2{x) > Ki{x) > 
Ki/2{x) which can be used for stringent numerical estimates for integrals containing these 
Bessel functions. 

We conclude this paragraph by noting that the x-space techniques allow one to compute 
sunrise-type diagrams in closed form in terms of elementary functions as long as one is dealing 
with odd-dimensional space-times. The resulting expressions are rather simple and can be 
directly used for applications. Having the complete formulas at hand, there is no need to 
expand in the parameters of the diagram such as masses or external momenta. 



5.2 Large powers of propagators 

In many physics applications one encounters large powers of propagators. Recent examples are 
the calculation of the corrections to — B^ mixing in perturbative QCD |182[I183] (see Fig. fT2|l 
or the large mass expansion for the contribution of charged scalars to the muon anomalous 
magnetic moment |184j . When calculating moments of the spectral density of such a diagram by 
using packages like MATAD for the automatic calculation of Feynman diagrams jl85[ I186[ I187j , 
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Figure 12: A non-factorizable diagram at NLO for — mixing 




Figure 13: Sunrise-type diagram with additional external momenta 

high powers of propagators are generated. Another example are sunrise-type diagrams with 
external lines at small momenta qi (see Fig. I13|). These diagrams appear in calculations when 
there are weak external fields as, for example, in sum rule applications |188| I189j or in special 
cases of high precision calculations in ChPT |107j . 

In momentum space larger powers of propagators are treated through recurrence relations 
based on the integration- by-parts technique |19()^ll9l]ll92j . In configuration space larger powers 
of the propagators give rise to larger values of the indices for the corresponding Bessel functions. 
We now describe a technique using recurrence relations which allows one to reduce large indices 
of Bessel functions |3n] ■ 

Note that the direct reduction of a sunrise-type diagram to a standard set of master integrals 
with the help of algebraic computer systems is rather time-consuming with present momentum 
space techniques. In practice the computation proceeds through the use of a table of integrals 
with given powers of the denominators. One would have to set up a three-dimensional table 
for a given total power N. The number of entries (even when accounting for the appropriate 
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symmetries) still grows as fast as A^'^ which is large for the large values of needed in some 
present applications. Within our method one first re-expresses the relevant integrals through 
a one-parameter set of integrals which are then solved explicitly. For large the number of 
entries increases as a first power of (the number of elements for the /o('?) basis is given by 
2[N/2] — 5 where [z\ is an integer part of z) which considerably reduces the time consumption 
in a computer evaluation. 

In order to show the applicability of our method we shall calculate three examples that 
were solved before with the help of momentum space techniques |136j . In Ref . |136j Broadhurst 
considered general three-loop bubbles B^. A subclass of these are the sunrise-type three- loop 
bubbles. In momentum space they read 

f d^k d^l d^p 
BN{0,0,ns,n,,n,,n,) = J ^3Z>(^i./2r(3 _ D/2))3 ^ 

m^"3 m^""* m^^s m^"^" 
X (182) 

((p + A;)2 + m2)"3 ((j9 + /)2 + ^2)n4 (^(^p_^^_^/)2_^^2)n5 (^p2_^^2)n6 ^ ^ 

with two propagators absent {rii = n2 = 0). Actually we choose the remaining indices 77,3, 
n4, ns, and ng such that the results become finite. We shall not only calculate the finite part 
but also the part porportional to e in order to be able to compare with |136j . Written in 
configuration space, the particular subset of bubble diagrams -Bat is given by 

B^iO,0,n„n,,n„n,) = y ( 1^5)^13^2 - m^^"3+n.+..+.a)-i2+6. ^ 

(x, m)D("*~i) (x, m)D("'-'-^) (x, m)D("«-^) (x, m)x^^+^dx. (183) 







It is clear that one will end up with integrals of products of four Bessel functions with non- 
integer indices and a non-integer power of x. We first discuss three examples in the next 
three subsections which will be followed by more general considerations on a general reduction 
procedure which allows one to reduce all integrals containing products of Bessel functions to 
the following set of two master integrals, 

U{r):= {Ko{0)Cd^ and L[{r) := {Ko{x))' C He'^'^m dx (184) 
Jo Jo 

where the index "Z" in L''^{r) reflects the logarithm appearing in the integrand. We will then 

add some considerations on the basic integrals L^lr) and L'-^{r). 

5.3 The example Bn{0, 0, 2, 2, 2, 2) 

We start with the example which is represented by the diagram in Fig. IT^ b). Each of the lines 
is modified (indicated by the dots on the lines) which means that instead of the propagators 
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Figure 14: The diagrams for 5^(0, 0, 2, 2, 2, 1), 5^(0, 0, 2, 2, 2, 2), and 5jv(0, 0, 2, 3, 3, 4). 



/^(x, m) we have to use 



One obtains = mx 



(27r)«(p2 + m2)2 2(27r)^+ 



5^(0,0,2,2,2,2) 



)l-2e 



(r(i + 5))3r(2-e) Jo 



(x/m)" 

— — — KAmx) 
2(27r)2-^ ^ 



We now can use the general formula 



'dK,{z) 
dv 



n-l 



= ±-n\ y 

2 ^ 

!/=±n k=0 



k—n 



Kkjz) 
k\{n — k) 



{0,1,...} 



;i85) 



(186) 



(187) 



to expand the Bessel function in a series with respect to its index. In case of K_^{z), however, 
the first derivative vanishes and we obtain K_^{z) = Kq{z) + 0(e2). Therefore, in expanding 



(r(l + £))3r(2 - e) ^ (l + (1 + 27E - 2 In 2 + 2 In 0^ + 0{e' 

= 2 (KoiO)' (l + (1 + 2 ln(e^-e/2)) e + 0{e' 



(188) 



we obtain 



5^(0,0,2,2,2,2) = 2{l + e) {Ko{0)' + (MO)' eHe^'^^md^ 

Jo Jo 



2{l + 6)L^{3)+AeL[{3). 



(189) 



We have compared our numerical result with the analytical result in |136j and have found 
agreement. In fact, from the transcendality structure of the results in |136j one can numerically 
obtain the appropriate coefficients multiplying the transcendentals. One finds |^ 



LJ3) = -— + — C(3), 
4v ; 32SV 

, 3 3/1 

ii(3) = j^ + jLU - 
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5.4 The example Bn{0, 0, 2, 2, 2, 1) 

In the diagram of Fig. IT^ a) one of the hnes is not modified. Therefore, we have to deal with 
one regular propagator factor 

D^^\x,m) = D{x,m) = i^^K,{mx) = i^^K,_^(mx). (191) 
In this case we obtain 

5^(0, 0, 2, 2, 2, 1) = (r(^^^))3r(2-e) I (^-^^^^ ^^1-^(0^'+'^^^- (192) 
Using Eq. p87|) one has 

Ki-eiO = MO - |^o(0 + 0{e') (193) 

and ^ g 

The result reads 

roo roo 

5^(0,0,2,2,2,1) = A{l + e) {Ko{OfK,{Oed^-'^e (MOf^d^ 

Jo Jo 

+8e {KoiOfK.iOeHe-'^^mdC (195) 
Jo 

This result is not yet written in terms of L^lr) and L'-^{r) and will be dealt with after introduction 
of the reduction procedure (The result for e = was given earlier in Sec. 13. 6p . 

5.5 The example Bn{0, 0, 2, 3, 3, 4) 

In order to demonstrate the power of the configuration space technique also in a more complex 
setting we finally choose the diagram in Fig. EJc) as an example. The modified propagators 
now read 

(x/m)^+^ 



D(-^\x,m) = (^M)'^V _i_,(mx), D^^\x,m) = , 

^ ' ' 8(27r)2-^ 1 v > ; 48(27r)2-^ 

Using the expansions 



K_2-s{'mx) (196) 



K^i^siO = K,{0 + -Ko{0+O{e'), K.2-e{0 = K^iO + jK.iO + ^MO+Oie') (197) 
^e obtain 

2-6-2e ^oo 

5^(0,0,2,3,3,4) = 3 K^^iOiK^i^eiO) K^2-e{0e^''d^ 



l + e 
192 



rco f roo 

/ K,{i){K,{i)fK,{i)Cdi+- KoiOiK.iOfK.iOed^ 
Jo 96 JO 



+ 96 70 ^°^^)^^^^^))^'^^+96yo (^o(^)) (^i(^))^'^^ 

+ -/ Ko{OiK,{0)K2m'He'^^/m- (198) 
96 Jo 
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5.6 The reduction procedure 

Especially in the last expression given by Eq. ()198j) there are many different integrals which 
differ from the basis L^lr) and L''^{r). However, the integrands can be reduced to integrands 
involving Kq[^) and Ki{^) only by using the relation 

11 — 1 

K40 = 2^i^„_i(0 + i^„-2(0- (199) 

After the first step, namely the expansion of Bessel functions for non-integer indices, the above 
relation establishes the second step in our reduction procedure. Finally, we use 

^KoiO = and 

= -liKoiO + m)) = -MO-^^m) (200) 
to perform the third and last step. For instance one has 

Jo Jo at, 

= -KoiO{Ko{Ofe]^ + fKo{0-[^{Ko{OfCd^ = 

= ^I^^KoiOiMOf^^ed^ + rfKoiOiMOfC-'d^ = 

= -3/ {Ko{OfK^{OCd^ + r {Ko{0)C~'d^ (201) 
JO Jo 

and therefore 

L^\r) = / {KoiOfmrd^ = - / {Kom'r'd^ = - (202) 

JO 4 JO 4 

The reduction formulas are given in general by 

Lt\r) = ^((r-m + l)L(r-i)(r-l)-(m-l)L(r~2)(^))^ (203) 



(r) = L^((r-m+l)Ll(-i)(r-l) + 4™"i)(r-l)-(m-l)L|('"-^)(r)) 

n — m + 1 ^ ' 



Tlim) 

and are coded in MATHEMATICA in order to automatically reduce to the master integrals 
After executing the second and third step of the recursion the results for the above examples 
read 

5^(0,0,2,2,2,2) = 2(l + £)L4(3)+4e4(3) + 0(£2), 
5^(0,0,2,2,2,1) = 2L4(l)+4£Li(l) + 0(£2), 

5^(0,0,2,3,3,4) = 1^4(3) -^L4(5)-^L4(7) + 0(e). (204) 
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Let us concentrate on the 0{e^) terms. By matching our numerical results to the results 
obtained by using the RECURSOR package 

5^(0,0,2,2,2,2) = -^ + lc(3) + 0(£), 
5Ar(0,0,2,2,2,l) = Ic(3) + 0(£), 

5^(0,0,2,3,3,4) = ^ + 0{e) (205) 

57o 

we can determine the master integrals Li{r) for a few values of its argument. One has 

^4(3) = + ^C(3), ^4(1) = ^C(3), 16L4(3) - 2L4(5) - ^4(7) = 1. (206) 

In particular the last relation is interesting because it implies a sum rule for the weighted 
product of three Bessel functions, namely 

/ KoiOiK.iO) K2m'dC=:^. (207) 
JO 3 

This surprising identity has been cross-checked by numerical integration. Note that for the 
0(e°) contribution only odd values of r appear as arguments in the master integral L4{r). 

5.7 Laplace- type asymptotics 



The two master integrals Eq. ()184j) are the basic set for the reduction procedure. A first 
estimate of the numerical magnitude of these integrals can be easily inferred from the asymptotic 
expansion of the integrals at large q, 



L[{q) = ^^|^(vl/(2g)+7i^-31n2)(^l-— ^ + 0(l/g2)) (208) 



where \E'(x) = r'(x)/r(x) is the logarithmic derivative of the F-function. The direct configu- 
ration space representation is the most convenient for numerical evaluation. To see this, we 
modify the last term of the integrals -^4(5) and L\{q) by introducing parameters k and ni, 
respectively. In doing this, we find that the relations 

= (^(2g) + - 3 ln2) (1 - (210) 
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with K = 0.97 {ki = 1.17) gives good results for the basis rational (resp. log-type) integrals with 
an accuracy better than 1% for all values of the arguments g > 1 (g > 3). For g > 3 (g > 5) 
the relative accuracy is better than 10^'^. The leading order asymptotic formulas for the basis 
log-type integrals are less precise for small g than its analogue in the rational case because 
the integrand in Eq. p84|) is not positive for log-type integrals. Obviously, an exact solution 
via recurrence relations in momentum space will be much more complicated than these simple 
asymptotic formulae. Because the cancellation of significant figures in the numerical evaluation 
of terms with different transcendental structure is a dominating feature of the momentum space 
recurrence procedure for large g, the use of asymptotic formulas is worth considering. 

5.8 Irreducible numerator: three-loop vacuum bubble case 

Loop integrals may have numerator factors which involve the loop momenta and cannot be 
expanded in terms of the denominator pole factors. In such a case one speaks of irreducible 
numerator factors. Momentum space techniques can run into problems when non-trivial nu- 
merator factors appear. Take for example the numerator factor {ki ■ ^2) for a n-loop bubble 
with n + 1 massive lines. For n < 3 the momentum space integral can be reduced to scalar 
integrals. For n = 3, however, the numerator factorintegral is no longer reducible |137j . The 
problem of irreducible numerator factors has a straightforward solution in configuration space 
by the integration-by-parts technique |141j . The starting expression involving the non-trivial 
numerator factor {ki ■ k2) (or any other scalar product of linear independent inner moments) 
corresponds to 

n*(0) = j D{x, mf {d^D{x, m)) (9^(D(x, m))) d^x (211) 

where the asterix denotes the fact that we are dealing with a non-scalar master integral. Careful 
use of the integration-by-parts identities |149| ll5Uj 

J d^,{D{x,m)---D{x,m))d'^x = 0, J D {D{x,m) ■ ■ ■ D{x,m)) d^x = (212) 

leads to the result 

J D{x,mf {d^D{x,m)) {d^'{D{x,m))) d^x = J D{x,mfDD{x,m)d^x. (213) 
For the last integral we have used (— □ + m?)D{x^m) = 6{x) and end up with 

J D{x, mYDD{x, m)d^x = J D{x, mfd^x - D{0, m)3. (214) 
The value of D{0,m) can be found for instance by integration in momentum space, 

U[0,m)- rn^+p2- (27r)^r(D/2) io p^ + ~ iAn)^/^'^' "^l^'- ^^'^^ 
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In the end we obtain 

n;(0) = -— n3(0) + -D(0,m)3 (216) 
where 113(0) is the scalar three- loop sunrise- type bubble diagram without a numerator factor. 

5.9 Irreducible numerator: four-loop vacuum bubble case 

Next we consider a four-loop diagram which appears as a second independent master integral 
V2 of the sunrise topology in the classification of Ref. jl37i|. In momentum space the second 
master integral V2 has the additional numerator factor {ki ■ k^)'^ as compared to the first scalar 
master integral Vi. The second master integral reads 

™ ... ^ f ih ■ k,fi27T)-^^d^k,d^k2d^k,d^k, 

^ J {ml + kl)im^, + {k2-k,)-'){mi + {ks-k2r){ml+{k,-ks)^){ml, + kiy ^ ^ 
Turning again to the equal mass case and using the configuration space representation, this 
integral can be written as 

n*(0) = J D{x,mf{df,d^D{x,m)){d^,d^D{x,m))d^x. (218) 

It is apparent that by using integration-by-parts techniques this integral cannot be reduced to 
scalar integrals and/or integrals containing d'Alembertians. The easiest way to evaluate such 
an integral is to compute the derivatives directly. This is done with the help of the relation 

(z-'Kiz)) = - (z-'-'K^+.iz)) (219) 

(cf. Eq. (IA12|) ). Eq. (|219j) can be iterated and gives results for arbitrary high order derivatives 
of Bessel functions Kx{z) in terms of the same class of Bessel functions with shifted indices and 
powers in z. For the first derivative we obtain 

o ^ m^^+^ Kx+i{mx) 

dMx,m) = -x,j^^^—^. (220) 

Since the resulting analytical expression for a given line of the diagram lies in the same class as 
the original line, the procedure of evaluating the integral is similar to the usual one. However, 
we cannot use the second derivative 

2 A-|-2 

d^d,Dix,m) = (2^)aTi(^^)a+i (^?/-^A+2(ma;) - ^Kx+i{mx)^ (221) 

directly under the integration sign. The reason is that the propagator has to be regarded as a 
distribution. There is a 5-function singularity at the origin which is not taken into account in 
Eq. ()22H) . Indeed, contracting the indices fi an u in Eq. 1)2211) one obtains 

d^df^Dix, m) = m^D{x, m) (222) 
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while the correct equation for the propagator reads (—9^ + m?)D{x,m) = 6{x). Thus, a 
straightforward evaluation of derivatives is valid only for x 7^ 0. The behaviour at the origin 
(x = 0) requires special consideration. In practice, to treat this case one should not use higher 
order derivatives but stay at the level of the first derivative. 

In order to deal with this case we introduce another master integral 

n**(0) = J D{x,m)d^D{x,m)d^D{x,m)d''D{x,m)d''D{x,m)d^x. (223) 

The relation between the two master integrals 114(0) and 114* (0) is found to be 

n:(0) = 3nr(0) - Im^fim - lm''D{0,m)\ (224) 
o o 

The quantity 114* (0) can be calculated with the explicit use of first order derivatives within our 
technique. The analytical result for the pole part reads 

A/4ll4W-m ^^^^3 4147205 ^^^V ^ ^ 

(A^ = ((47r)^"'^/r(l = e)Y, cf. Eq.dSn}) and the 5-expansion in the form 

7V'4n**(0) = m^°( - 0.375^-^ - 1.923611e"^ - 5.474103£"2 - 11.90356£"^ 

-27.99303 - 104.53845 - 663.61235^ - 3703.241^^ + 0{e^)) . (226) 



Since the results for 114(0), Il4*(0), and 1^(0, m)^ are known, Eq. ()224|) can be used to obtain 
the final result for the original integral. 



Af^UliO) = m^°(- 1.6875£"^-7.8125£"^- 21.20964£"=^- 44.769555"^ 

-97.07652 - 290.92345 - 1719.8095^ - 8934.731^=^ + 0(£^)) (227) 



which again verifies the result given in Ref. |137j . 

Differentiation of the massive propagator leads to expressions of a similar functional form 
which makes the configuration space technique a universal tool for calculating any master 
integral of the sunrise topology. This technique is also useful for finding master integrals. 
Indeed, new master integrals appear when there is a possibility to add new derivatives into the 
integrands which cannot be eventually removed by using the equations of motion or integration- 
by-parts recurrence relations. But once again: without explicit inclusion of the 5-function only 
one derivative is allowed. Otherwise one misses tadpole parts of the result. Therefore, the new 
master integral should contain just one derivative for each line except for one line. This allows 
one to enumerate the number of non-scalar master integrals, i.e. master integrals including 
non-trivial numerator factors. For instance, in the five-loop case (six propagators) there will 
be only two non-scalar master integrals. 
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Figure 15: spectacle topology diagram in two different representations (a) (left hand side) and 
(b) (right hand side). The configurations space points 0, x and y are indicated. 

6 Generalization: Still within reach 

Apart from the diagrams with sunrise-type topology treated in the previous sections there are 
more involved topologies close to the sunset-type topology that still allow for quite a simple 
numerical evaluation. We conclude this report by discussing several such cases which are still 
in reach of the methods presented here. 

6.1 Generalization to the spectacle topology 

In this subsection we give a formula for a more general topology which occurs when one prop- 
agator is removed from an initial three-loop, tetrahedron-shaped bubble diagram. In the origi- 
nal classification of Ref. |in8j these diagrams are class E diagrams belonging to the spectacle 
topology (see Fig. The formulas obtained in this subsection are efficient for the numerical 
integration of diagrams of the spectacle topology. We have, however, not been able to find an 
analytical solution to the problem. The main obstacle of generalizing the configuration space 
technique to a general multi-loop diagram is the angular integration. The configuration space 
technique proved to be rather successful for general diagrams in the massless case but it 
brings no essential simplification in the general massive case (see e.g. Ref. [2S])- However, for 
special configurations the angular integration can be explicitly performed with a reasonably 
simple integrand left for the radial integration. Diagrams of the spectacle topology such as the 
ones shown in Fig. ^1 are examples of such a configuration. 

The configuration space expression of a spectacle topology diagram can be written in a form 
suitable for our purpose (see Fig. IT5la)). 



The key relation which leads to the simplification of the configuration space integral with 
spectacle topology is an addition theorem for Bessel functions in Eq. ()A17|1 . allowing to perform 
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the integration over the relative angle in the propagator D{x — y,m). Using Eq. ()A17|) for 
Z = and substituting m = e*'^/^, we can perform an analytic continuation in order to 
obtain a relation for the modified Bessel functions K and /, 

= 2-r(A)f;(A + i):^%Mc,^(„Bv,) (229) 

A:=0 P ^ 

where we have R = -^/r^ + — 2rpcos(p, 

I^{z) = e-"^'/VA(e"'/22) for - vr < arg^ < ^ 



2' 



I^(z) = e^^'^'/^Jxie-^'^'^^z) for ^ < arg^ < tt. (230) 



Using the orthogonality relations for Gegenbauer polynomials in Eq. ()B3|) . the sum disappears 
after integration over the relative angle and only one term contributes. We obtain 

J = 2 T{X)}_^{X + k)—-^ -J — J C,,{cosif)dQ^ 

^ 2^r,.,A%Mi_?^^c^„ 

- r>,, (231) 

where the first equality is a consequence of the orthogonality relation with the trivial factor 
Cq (1) = 1. This result allows one to write down an expression for any spectacle- type diagram 
in the form of a two-fold integral with a simple integration measure, 



Jo Jo 



X 



POO POO 

/ D{x,myix{x)x^+^dx / D{y,m)^Kx{y)y^^^dy (232) 

Jo Jx 



where 9{x) is the standard step-function distribution 

Note that the integration measure D{x,myx'^^~^^dx allows one to perform the integration 
by using efficient integration routines for numerical evaluation. The form of the weight function 
is close to e~"'^x" which suggests the use of Laguerre polynomials 

Lnix) = ^— {x-e-^) (233) 

as a convenient choice within the Gaussian numerical integration method. In fact, any modified 
propagator (with any power of the denominator) can be used as a factor in the integration 
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measure D{x,mYx'^^^^dx. This fact makes the representation universal and useful for the case 
of higher powers of denominators of the lines associated with pairs (x, 0) and (y, 0) of space- 
time points. If the angular structure of the diagram is preserved, the generalization to higher 
loops in the expressions for the radial measures is straightforward. 

6.2 Example for the spectacle diagram 

To show how this technique works, we give an example of the evaluation of a spectacle diagram. 
Consider an integer dimension space-time which, without loss of generality, can be chosen to be 
two-dimensional (an odd number of dimensions is trivial because, as we have seen earlier, the 
propagators degenerate to simple exponentials). The spectacle-type three-loop diagram can be 
obtained in a closed form. Indeed, in momentum space representation we have 

for the basic spectacle diagram S with 11 (p) and the mass M of the connecting propagator kept 
different. After the substitutions p = 2msinh(^/2) and t = we obtain 

where 7 = 1 — M"^ /2m?. Performing the integration, we finally obtain 

S{M) = f^^'^-f^^^^ (236) 



with ti^2 = 7 ± -\/7^ — 1 and 



_ 3tL.3(lA) - (i + 7)C(3) 



Lis (2;) is the trilogarithm function 



Li3(^) = E 4 kl < 1- (238) 
k=i 



For M = 2m (i.e. 7 = —1) the integral in Eq. ()235|) simplifies and one finds a simple answer in 
terms of the (for the present context) standard transcendental numbers In 2 and C(3), 

^(2m) = -^fIc(3)-ln2V (239) 



Anm^ \8 

For the equal mass case M = m we obtain a result in terms of Clausen's trilogarithm Cl3(27r/3). 
As one can conclude from these results, the conjugate pair of the sixth order roots of unity. 
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exp(±27ri/3) play an important role in this case again in accordance with the general analysis 
in The origin of the appearance of the sixth order roots of unity as the parameters of 

the analytical expressions of the diagrams lies in the mismatch of masses along the lines of the 
diagrams. However, the exceptional case M = 2m, where one line has twice the mass of the 
other lines (which results in the drastic simplification), also keeps us within the set of the sixth 
order roots of unity. The key parameter in this case is simply the natural number 1 which is 
definitely one of the sixth order roots of unity. 

Turning to the configuration space representation in Eq. ()232j) we find 



S{M) = / xKo{mxydx / yKQ{myYdy x 
JO Jo 

(Ko{Mx)Io{My)e{x - y) + K,{My)h{Mx)e{y - x)\ . (240) 



for the basic spectacle diagram. An explicit numerical integration of Eq. ()240|) shows coincidence 
with the analytical result in Eq. ()236|1 which we have checked numerically for arbitrary values 
of M. In this example the analytical result has a rather simple form. This is not true if 
high powers of the denominators enter. Then the corresponding one-loop insertions are rather 
cumbersome and an explicit integration in configuration space is more convenient. 

6.3 About the occurrence of Clausen's dilogarithms 

We want to have a closer look at how the square of Clausen's dilogarithm Cl2(vr/3)^ can emerge 
at the level of spectacle topology diagrams. As transcendental number, Clausen's dilogarithm 
Cl2(vr/3) characterizes the analytical results for three- loop bubbles. Its presence was discovered 
in the impressive treatise of David Broadhurst on the role of the sixth order roots of unity for 
the transcendental structure of results for Feynman diagrams in quantum field theory 

We consider the spectacle diagram in the form shown in Fig. [T^T b). Here the expression for 
the generalized middle line is a product of the one-loop propagator and the standard particle 
propagator. We express the one-loop propagator using a dispersion representation with the 
spectral density p(s) and obtain 



+ + Jim? s + p^ Jim^ s — \p'^ + s + p^ 

We take only the first term which is sufficient for obtaining the result we are aiming for. In 
this case the integral becomes independent of p^ and can be considered separately. One has 

1 r Pis)ds ^242) 



pz _|_ jy^z J s — m? 
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which leads to the sunrise diagram after the two other hne shown in Fig. IT^T b) have been added 
with a normahzation factor given by the integraL One factor Cl2(7r/3) results from integrating 
the overall sunrise diagram which is composed of the propagator {p^ + m^)~^ from Eq. ()242|) . 
The two other lines of the diagram are shown in Fig. lTH b). The second factor Cl2(vr/3) has to be 
found in the normalization factor given by the integral in Eq. ()242|) . Note that the very structure 
of this contribution - the square of a number which first appeared at the lower loop level - 
suggests a hint for its search. It should emerge as an iteration of a lower order contribution in 
accordance with the iterative structure of the R-operation which provides a general framework 
for the analysis of multi-loop diagrams. The following consideration confirms this conjecture. 
Consider the quantity 

N=r (243) 

and take p{s) to be the spectral density in D-dimensional space-time |SHI, 

(s-4m2)^-i/2 

P(^) = 24A.i,..V2r(A + l/2)V^ ' ^ > 2m. (244) 

Next consider the first order contribution of the expansion in e near the space-time dimension 
Dq = 2. The expansion in X = e results in 

(s-4m2)^-i/2 1 / fs-Am^\ 

Therefore, the relevant first order term in e reads 

A.p(.) = H^^'^-^')/-^') (246) 
2ttJs{s - Am?) 

where /i = m has been chosen for convenience. Next we change variables according to 

= 2m cosh(^/2), t = e'^ (247) 

to obtain 

. / 9 9/ , XX (-Int + 21n(l -t))t 

A,p(4m^ cosh^(e/2)) = ^ 3^^2(1 1, 2) • (248) 



For the quantity in question we find 



/4m2 S — m? 271711? Jo 1 + t + 

The roots of the denominator of the integrand in Eq. ()249|) are now given by 4 = exp(±27rz/3) 
which again is a conjugate pair of the sixth order roots of unity. After integrating this equation 
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we readily find 



^'^ = -^^(-f'"3 + Im(Li.(e-/.)_Li,(e--/3 



27rm2V3 V 3 

^ ^ClJ^U^lnS). (250) 



Using the relation 
we finally obtain 



irm^VS V V 3 ; 6 

Therefore, in the first order of the e expansion of the spectacle diagram we indeed found the 
remarkable contribution proportional to Cl2(7r/3)2. In our calculation it emerges naturally as 
the iteration of the lower order term (42j . Originally, the presence of this contribution has been 
conjectured and confirmed in |llUj by direct numerical computation of the finite part of the 
general three-loop bubble in four- dimensional space-time. 



6.4 The insertion of massless irreducible loops 

One can consider a sunrise-type diagram with one or more irreducible loops by which we mean 
a generalized line more complicated than an ordinary standard propagator or a (large) power of 
them. Consider, for instance, the replacement of a line by a subdiagram of fish-type topology 
as shown in Fig. ^| The configuration space technique leads to a numerical solution in this 
case because we can replace one of the propagator factors D{x, m) by the two-point propagator 
in the irreducible subdiagram, using the fact that 

n(x) = j D{x, ^/s)p{s)ds (253) 

where p{s) is the spectral density of the subdiagram. If the subdiagram is a massless fish-type 
diagram, the spectral density is given by 

p{s) = 6C(3)5(s). (254) 

The known result in momentum space can be restored by using the dispersion relation (j^OJ. 
On the other hand, we obtain 

D{x) = 6C(3) J D{x, ^)5{s)ds = 6C(3)L>(x, 0) (255) 

which up to a factor is the usual massless propagator. 
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Figure 16: Irreducible loop replacing one of the propagators 



6.5 The insertion of massive irreducible loops 

Insertions of massive irreducible loops were used analytically for the calculation of corrections to 
baryon correlators jMl IHEl ISHj and numerically for the mixing of heavy neutral mesons [H] . The 
massive insertion can be incorporated mainly because the two-loop diagrams are well analyzed 
with any mass configurations |198| I194| I195| I196j . For the semi-massive fish-type one-loop 
diagram with one massive and one massless line we can use 

But even in the case with two different masses a result for the spectral density is available, 

pis) = ^^J^i^ (^hixi) + Li2(x2) - ^ In(-xi) ln(l - xi) - ^ ln(-X2) ln(l - ^2)^ (257) 
where 

2m2 

xi = , 

ml + ml — s + J (mf + m2 — s)^ — Amlml 

2ml 

X2 = . ' (258) 

mf + ml ~ s + J (mf + — s)^ — Am^m^ 
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7 Summary and conclusion 



To conclude, we have presented a review of configuration space techniques for the calculation 
of sunrise-type diagrams. We have shown that the singular or pole parts of any sunrise-type 
diagram are calculable analytically in the simplest possible manner. For the finite parts we ob- 
tained one- dimensional integrals of well-known functions which is very convenient for numerical 
evaluation. We have presented many sample calculations of multi-loop sunrise-type diagrams 
and have compared them with momentum space results in the literature when available. We 
have found agreement in every case. The agreement provides for a mutual check of the re- 
sults which have been derived using very different methods. In the asymptotic analysis have 
dealed with different expansions. We can conclude that any kind of expansions is given as an 
expansion in parameters of the propagators and thus leads to the expansion coefficients in a 
straightforward way. Finally, extensions to non-standard propagators and other exotic settings 
figured out to be feasible as well as some considerations for slightly different topologies. All in 
all, we can stress again that the benefit of configuration space techiques is the fact that there 
is 'almost' no integration in sunrise- type diagrams, as Eq. ^ shows. 
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A Properties of Bessel functions 



Since Bessel functions play a crucial role in our calculations we collect a number of definitions 
and properties of Bessel functions, mostly taken from the handbooks of Watson (1944) |133j . 
Prudnikov, Brychkov and Marichev (1990) |134j . and Gradshteyn and Ryshik (1994) 

Bessel functions and Hankel functions 

Ordinary Bessel functions are solutions of the Bessel differential equation 

^,fZ^ + + {z' - u')Z,{z) = (Al) 

az"^ az 

[u need not be an integer). There are two classes of solutions Z,y{z). Bessel functions of the 
first kind Ju{z) are nonsingular at the origin 2 = while Bessel functions of the second kind 
Yy{z) are singular at 2; = 0. Combining Bessel functions of both classes, we end up with Hankel 
functions 

H^{z) = J,{z)±iY,{z) (A2) 

(instead of Hf{z) the notation H^^^ = H+ and = R- are quite common). 
Bessel functions of the first kind have the series expansion 

Including a factor occuring in practical applications, we can also write 

2) ^''^'^-l',k\r{u + k + i)-r{i + u)[^ 1 + Z.+--J- ^^^^ 

T{z) is Euler's Gamma function. Functional equations are of help in order to relate Bessel 
functions of different degree. They are given by 

zZ^^i{z) + zZy+i{z) = 2iyZ^{z), Z^^i{z) - Z^+i{z) = 2-^Z^{z) (A5) 

where Z is any of the Bessel functions J, Y, or H^. As a consequence we obtain 



dz^ 



{z-'^Uz)) = -2-V,+fc(z). (A6) 
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Modified Bessel functions 



We can continue Bessel functions into the complex plane, ending up with modified Bessel func- 
tions of the first kind Iu{z) and modified Bessel functions of the second kind Ky{z)^ sometimes 
also known as McDonald functions. Iv{z) and Kjj{z) are solutions of the modified Bessel dif- 
ferential equation 

(AT) 



dz^ dz 
The modified Bessel function of the second kind 



dZyi^z) 2 2\ry / \ n 

- [z - U )Z^{z) = 



KJz) 



■Kl^y{z) - I^{z) 



Sm TTZ/ 



(A8) 



2 sin(7rz/) 

can be expressed by the modified Bessel function of the first kind Iu{z) with series expansion 



V CO 



,5 /c!r(i. + A; + 1) 



(A9) 



Therefore, one has 



:rl KAz) 



\-v\2 



r(i - v) (z 



r(i + z/) V2 



2v 



+ 0{z\z 



4 ^2+2v\ 



(AlO) 



For Eq. ()A5|) one can we derive functional equations for the modified Bessel functions, 

zlu-.i{z) - zlu+i{z) = 2ul^{z), h-iiz) + h+i{z) = 2-^I^{z), 

zK,^,{z) - zK,+^{z) = -2pK,{z), K,.^{z) + K,^i{z) = -2-^K,{z), (All) 

dz 

For the differentiation of the modified Bessel function including an appropriate factor we derive 



d 



-{z-Mz) 



-z-'K.+Az). 



(A12) 



Bessel functions can also be differentiated with respect to their index, if this index is given by 
an integer. For Ky{x) we for instance obtain 



dK,{z) 
du 



1 



n-l 



k=0 



k\{n — k) ' 



ne {0,1,...}. 



(A13) 



Finally, the modified Bessel function Ki, and the Hankel function are related by 

K.{z) = ye--/2if,+ (zz). 



(A14) 
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Asymptotic behaviour 

Asymptotic expansions determine the behaviour of the modified Bessel functions of the second 
kind and the Hankel functions H~ when the arguments become large. For the two functions 
just related to each other we obtain 



TT 



2z 



1/2 



N-l 



and 



\7TZ, 



N-l 



r(z/ + n - 1/2) 
n!r(z/ - n - 1/2) 

-l)^(j.,iV)^ 



(2z) 



TV 



(A15) 



(A16) 



where 9 G [0, 1] is chosen appropriately. For the asymptotic behaviour of Iu{z) see Eq. ()C16|) . 
Addition theorem 

The addition theorem for Bessel functions is given by 

= 2^m-^r(.) £(. + fc):^l±4;^^^^Cncos^) (A17) 

where are the Gegenbauer polynomials (cf. Appendix El), Z is any of the Bessel functions 
J, y, or H^, R = i/r^ + p2 _ 2rpcos(p, and r > p. For r < p the arguments of the Bessel 
functions on the right hand side of Eq. ()A17jl should be interchanged. 



B Properties of Gegenbauer polynomials 

Gegenbauer polynomials can be generated using the characteristic polynomial 

oo 

{t^ - 2tx + = 5]rc;;(a;). 

n=0 

In particular, we have Cq{x) = 1, C'({x) = 2vx and 

{n + l)C„%i(x) = 2(n + v)xC:{x) -{n + 2v- l)C:_,{x). 
Gegenbauer polynomials obey the orthogonality relations |lllj 



C^{Xi ■ X2)CI[{X2 ■ Xs)d^l2 



C^iXi ■ X3) 



2tx 



(Bl) 



(B2) 



(B3) 



r(z/ + i)n + z/ J r(z/ + i) 

where Xi are unit four-vectors and dQi is the angular part of the integration measure. Finally, 
we have 



c:(i) 



r(n + 2z/) 
n\T{2u) 



(B4) 
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C Cuts and discontinuities 



Functions which are continuous at least on parts of the real axis can be continued to the whole 
complex plane. However, if singularities occur on the real axis or elsewhere in the complex 
plane, this continuation need no longer be unique. Examples which occur in the context of our 
calculations are logarithms and polylogarithms but also powers with non-integer exponent and 
Bessel functions. We will deal with these special cases in this appendix. 

The case of logarithms 

For X — the function \n{x) has a non-removable singularity while for negative real values 
it is not defined. However, we can continue to the negative real axis by following a path 
which circumvents the origin. It makes a difference whether we turn around the origin in the 
mathematically positive or negative sense. If we for instance start at the real value x and turn 
around in the positive sense to continue to a value at —x, we can parametrize this by writing 
z — xe"^ where (p starts at (p — and increases until it reaches the negative real axis for (p — tt. 
By doing so, we obtain 

ln(xe+^'^) =lnx + i7r. (CI) 

However, we can reach the negative real axis as well by starting with (p again but going to 
negative values, reaching the negative real axis for (p — — tt. In this case we obtain 

\n{xe-'^) =lnx -tTT. (C2) 

Obviously, the value "ln(— x)" is not unique. Indeed, we can even think of the continuation to 
the complex plane as a sheet which winds up and rises each time we turn around the origin, 
the so-called Riemannian sheet. Nevertheless, for only one turn we can calculate the difference 
between the values. This quantity is known as discontinuity, and it can be calculated actually 
for every complex value z except for singular points. The definition is given by 

Disc/(^) = /(^e+*) -/(^e-*), e±^° = lim e±^^ e > (C3) 

where z = is the singular point of the function f{z). If we apply this definition to the case of 
a logarithm with negative argument, we obtain 

Disc In(-x) = ln(-xe+^°) - \n{-xe'''^) = ln{xe-''') - \n{xe+''') = -2me{x) (C4) 

where 9{x) is the step function {6{x) = 1 for x > 0, 6'(a;) = for x < 0). Note that we have 
been careful in replacing — a;e^*° by xe^*'^. The discontinuity of powers of logarithms can be 
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calculated as well, for instance 

Discln^(-x) = ln^(xe-^^) -ln^(xe+^'^) 

= ((Inx-iTr)^ - (lna; + i7r)^)6'(a;) = -Amlnx e{x). (C5) 
For later convenience we divide the discontinuity by the factor 27ri to obtain 

-^Disc In(-x) = -0{x), ^Disc In^(-x) = -21nx 0{x), . . . (C6) 

The case of poly logarithms 

Polylogarithms can be defined iteratively by 

UJx) = r ^^^^^=4^dx', UAx) = - ln(l - x). (C7) 
Jo X 

For all the polylogarithms a discontinuity occurs for x > 1. We obtain 

Lii(xe±^°) = - ln(l - xe^'°) = - ln(l + xe^'"") = (- ln(x - 1) ± in) e{x - 1) (C8) 

where we have extracted a factor e^'^ from the argument, and therefore 

-^Disc Lii(a;) = e{x - 1). (C9) 

In order to calculate the discontinuity of Li2(x), for x > 1 we split the integral into two parts, 

1 ln(l -x') , , r ln(l -x') , , 



Jo x' Jl 



x' 



^ u,(i)- r'^^^d.i.ir^. (CIO) 

Jl X Jl X 

The real parts cancel out when taking the difi^erence needed for the determination of the dis- 
continuity. Contrary to this the imaginary parts add up and one obtains 

1 r dx'. ^ . , /-^ dx' 



1 \ r dx dx 

DiscLi2(x) = / — DiscLii(x') = / — =lnx^(a; — 1), in the same way 

2m 27ri J x' Ji x' 

-^DiscLiafx) = ^— / ^DiscLi2(x') = / ^ Inx' = In^ a; ^(x — 1), in general 

ZTTZ ZTTZ J X' Jl X' 2 



1 [ dx' ^ . , ,^ /-^ dx' 
-^DiscLi„(a;) = In""^ x ^(x - 1). (Cll) 
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The case of non-integer powers 

For powers with exponents close to integer values such as {b — ax)~^ the calculation of the 
discontinuity is closely related to the calculation of the discontinuity of the logarithm because 

= l + e\nx + 0{e^). (C12) 

Similarly, the discontinuity of (6— aa;)~" can be obtained in the more general case in closed form. 
In order to calculate the discontinuity of a power (6 — ax)~" with non-integer exponent —n (we 
include the minus sign for later convenience) we first have to ponder where the discontinuities 
can appear. Obviously, this is the case for 6 — ax < 0. Taking a > and 6 > which is the 
case for all applications in this report, we obtain x > h/a. Therefore, we obtain 

Disc (6 -ax)-" = (6-ax+'°)-"-(6-ax"'°)-" = (6 + axe"*")-" - (6 + axe+'")-" 

= (e+^™(ax-6)-"-e-^™(ax-6)-")^(ax-6) = 2i sin(7rn)(ax - 6)-"^(ax - 6). (C13) 

Using 

sin(7™) = y = - sin(-7m) (C14) 

i [n)\. (1 — n) 

where r(x) is Euler's Gamma function, we finally obtain 

1 T^- /u \-n {ax -hy6{ax -h) 

Disc ih -ax) = , ' : — - (C15) 

2m ^ ' T{n)T{l-n) ^ ' 

Note that for negative integer values of n the discontinuity vanishes. 
The case of Bessel functions 

The final item about ambiguities found in the continuation into the complex plane we want to 
mention here is the so-called Stokes' phenomenon which occurs for the asymptotic expansion 
of Bessel functions and is mentioned already in the handbook of Watson |133j . The asymptotic 
expansion of the modified Bessel function of the first kind reads 

^ - (-1)^ r(z/ + fc + l/2) 

^ v^to (2^)' r(z/-A;+l/2)A;! 

-.±(.+i/2Ki oo^ 1 r(z/ + A; + l/2) 

where the plus sign is valid for — 7r/2 < argz < 37r/2 while the minus sign has to be taken for 
— 37r/2 < argz < 7r/2. The surprising fact is that there are two different asymptotic expansions 
available at every point of the complex plane, except for the arguments ±i7r/2. Therefore, in 
continuing Jj,(±ix) = lyixe^™!'^') to the complex plane, the choice is still unique. 
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D Singular contributions up to four loops 



In the following we present complete results for the singular parts of sunrise-type diagrams 
with arbitrary masses up to four-loop order. The results are given in the MS-scheme in the 
Euclidean domain. 

U{{p, nil, 1112) = -, 
e 

UKp, mi, m2, 7713, m^) = ^ E ^^i^^'j + E "^i ~ E "^i + 2 E "^i "^i - 3^0 ) + 



6£3^ ' ^ 12£2 . 

I + V E ^?(3 - 2^0 - 9 E <(5 - 2^0 + 6 E m>2(20 - 24^^ + + e^^.,, , , 
72£ \ , , 



1 

n^(p,mi,m2,m3,m4,m5) = ~24£4 ^ mlm]ml + 

E ~ E(™^™? + ^) + 2 E ^l^^j^U^ - 

'2p^ E "^i - 6/ E + 2 E + 3p' E - 



1 

'288£2 



-6 E("^i + - 4^i) + 12 E m^im^jml{15 - 2Ui + 2^' + Uilj) 

^ 3/ + 2p^ ^ m2(35 - 24£i) - 18^^ ^ m^(21 - 8^^) + 2 E "'l'(77 - 24^^) 



1728e 

+V E mlm]{n - SUi + 8£2 + 2AUlj) - 216 E(^i ^? - mlrri^)U 
-18 E(^i^? + ^i^J)(49 - 56£i + 4£2 + I2£i£^) 

+24 m^im^jml{105 - ISOii + 30 + QO^i^j - 2£^ - ISi^ej - I2iii/k) J (Dl) 

i¥=j^k J 

where ii = ln{mf / fi'^) . The indices i, j, and k run over all mass indices. One can check that 
the general results listed in this Appendix reproduce the equal mass results in the main text. 
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E Subtraction terms for the small x singularities 



The leading singularity at small x is given by the massless propagator of the form 

The next order of the small x-expansion for the propagator D{x,m) is explicitly given by 



This term is suppressed relative to the first term by one power of at small x in four- 
dimensional space-time (however, this is not the case for two-dimensional space-time with 
A = 0). The term 

is further suppressed by one power of a;^ at small x. Therefore, the full subtraction of the 
three terms gives a rather smooth behaviour at small x which is sufficient to obtain a regular 
integrand for the numerical integration. 



F Analytical results for the four-loop sunrise diagram 

In this appendix we present some more details of our calculations for the four-loop sunrise 
diagram. For the analytical evaluation we take the last two terms of the integrand from Eq. (fTSj) . 
One has to integrate a product of two Bessel functions multiplied with powers of x which can 
be done analytically. The explicit expression for the e-expansion of that part of the integral 
which is evaluated analytically reads 



nf^(O) = 



5 35 4565 58345 



2e^ 3^3 144^2 8QAe 
1456940638037 1 709971^ 38 577r4 ^ 2525968C(3) 



7779240000 24 10 105 

/ 55171475321621447 245750977^ 129 2 5 377r^ 

^\ 1633640400000 ^ 144 175 

6752474831C(3) . , , , A 

+ ^ + 1653077^^ 3 + 59508C 5 ]e 

44100 ' 7 

/ 10610679621089130529 9278194977^ 42901137597r^ 2259177"^ 
+ + 



68612896800000 864 110250 14 
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/ 5963907632629558995931 13252040337r2 
^\ 14408708328000000 ^ 5184 
4643790856997r^ 4852923l7r6 



6615000 2205 
312138383154103C(3) 72850437r2C(3) 4 238229084C(3)2 



+ + 435297rX(3) - 



277830000 6 ' 105 

+ ^-^^^^^^^ + 2479507r\(5) + 1190160C(7)).^ + 0(.)^|. (Fl) 

This expression shows the real complexity of the calculation which reveals itself in the structure 
of the results. The main feature is that the terms cannot be simultaneously simplified to all 
orders in e. By a special choice of the normalization factor one can make the leading term 
and, in fact, even all pole terms simple, but then the higher order terms contain rather lengthy 
combinations of transcendental numbers that are not reducible in terms of standard quantities 
such as the Riemann (^-functions. Note also that the rational coefficients of transcendental 
numbers are very big and there is a huge numerical cancellation between the rational and 
transcendental parts of the answer (see also the discussion in Ref. pIU]). 



G Integrand for the numerical integration 

For the numerical evaluation we take the first three terms of the integrand from Eq. (fTSj) . One 
has to integrate them numerically as there is a product of three or more Bessel functions which 
is too complicated to be done analytically. To find the e-expansion of the integral one has 
to first expand the integrand in e. The expression for the e-expansion is quite lengthy. We 
therefore give explicit results only for e = 0. For this part the integrand for the numerical 
integration over z = mx reads 

™„m/ X fi/ , ,2 ,s 384 384 768/ 24 480/ 576/^ 
nf X = 66 - 108/ - 144/2 + 192/3 + — ^+ + ^ + 

y z^ z^ z^ z'^ 

111^2 ^47^ 2 ^g5 4 

+ mi^z^ + 2Afz^ + 24/^2 ' 



16 2 32 

201/2^ 9/2^4 75^6 405/;z6 531/2^6 

24/^2^ + 12/*^^^ + 



16 2 512 128 64 

15/3^6 9/4^6 375^8 ^2^l^& 315/2^8 5^^3^8 3^4^8 
+ + -^TTTT - + —737^ + 



2 4 2048 1024 256 64 16 

1875^1° 375/^^0 225/ V 15/ 3/V0\ , , 

H \ \ \KAz) 

131072 8192 4096 512 512 / ^ ' 
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./-512 384 768/ 24 288/ 384/^ , 

— — + — + — + 52z + 120/z - 64/^2 

\ z° z-^ z-^ z z z 

_ 21/^3 + 48/2^3 _ 24/3^3 + !^ _ 1^ + 9/2^5 _ 3;3^5 

8 32 16 



125^^ 75lz^ 15/V /^^^\,^,,2 A28KAzf 

\ ]Ki(zY + P--. 

512 128 32 8 / ^ ^ z^ 



Here / = ln{e1;z/2), z = mx, and je = — r'(l) = 0.577. . . is Euler's constant. As shown in 
Fig. 13 the plot of this function as well as the shapes of the corresponding functions in higher 
orders of e are very smooth and quite similar. The analytical expressions for higher orders in 
e, however, become much longer. Note that the new functions fn{z) first appear at order e^. 

The smoothness of the zeroth order integrand as shown in Eq. (jdljl implies that the numer- 
ical integration is quite easy to execute. Because the integrand vanishes exponentially for large 
values of z and has no singularities of the kind 2; In 2; for small values of z, the integration can in 
principle range from to 00. However, for practical reasons we had to instruct Wolfram's Math- 
ematica system for symbolic manipulations (which we used for all of our calculations presented 
here) that the integrand vanishes for z = 0. On the other hand, the asymptotic expansion of 
the integrand together with the integration measure is dominated by the part 

^^z''ln\e^z/2)K,{z)z'dz {K,{z) ^ ^e'^ for z ^ 00). (G2) 

Integrated from A to 00, this part gives a contribution 

67r^m^ 



A25/2in4(e^A/2)e-^ (G3) 



Hf -"(0) = 27T^ / Hf '"(x)x=^c/x ^ 27c^ / UT'^{x)x'^dx + A^''/Hn\e^A/2)e'^ (G4) 

Jo Jo 512 



512 

and terms which are of subleading order. Therefore, the result can be well estimated by 

512 

and A can be adjusted to any desired precision. 

A possibility to avoid any cutoff is to change the integration variable such that the interval 
[0, 00] is mapped onto [0, 1]. Then the integration can be done numerically with the additional 
information that the integrand vanishes identically at both end points. Possible transformations 
of this kind are for instance z = ln(l/t) or z = {1 — t) /t for t G [0, 1]. 
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